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Introduction 


The  analysis  and  understanding  of  short-period  regional 
phases  from  explosions  and  earthquakes  is  a  major  theme  in  CTBT 
issues.  If  a  total  ban  on  underground  testing  is  to  be  ade¬ 
quately  monitored,  then  the  means  must  be  available  to  discrimi¬ 
nate  between  explosions  and  earthquakes  down  to  very  small  seis¬ 
mic  magnitudes  and  at  relatively  short  distances.  Small  seismic 
magnitudes  ensures  that  short-period  data  must  be  examined  and 
that  high  frequency  wave  propagation  within  various  crustal 
structures  must  be  understood  to  unravel  the  source  effect  from 
wave  propagation  effects. 

A  major  effort  in  this  research  has  been  in  the  computation 
<  short-period  synthetic  seismograms  to  explain  data  from  local 
^nd  regional  earthquakes.  In  particular,  Section  1  describes  a 
study  of  very  shallow  foreshocks  and  aftershocks  associated  with 
the  1968  Meckering,  Australia,  earthquake.  Vertical  component 
short-period  seismograms  were  modeled  in  order  to  investigate  the 
excitation  of  depth  phases  within  the  P  arrival  and  the  excita¬ 
tion  of  short-period  Rg. 

The  purpose  of  this  study  was  to  find  accurate  depth  discri¬ 
minants  in  the  local  seismograms  themselves  rather  than  relying 
on  standard  earthquake  locations.  This  will  be  important  in  ana¬ 
lyzing  data  from  a  sparse  national  network  to  determine  source 
type  and  depth.  Depth  is  an  important  parameter  in  source  dis¬ 
crimination  since  the  likelihood  that  an  event  is  a  presumed 
explosion  is  small  if  the  event  can  be  accurately  located  below 
several  kilometers.  In  the  Meckering  study,  both  the  existence 
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of  sP  in  the  P  wave  form  as  well  as  large  Rg/S  ratios  served  to 
produce  very  accurate  depth  estimates  from  single  station  data  at 
ranges  of  60  to  95  km  from  the  source.  The  results  of  this  study 
yielded  several  interesting  conclusions  concerning  the  nature  of 
wave  propagation  in  shield  areas  and  in  the  tectonic  implications 
of  extremely  shallow  faulting  during  the  Meckering  event.  Ancil¬ 
lary  work  was  performed  on  the  teleseismic  and  static  uplift  data 
associated  with  the  Meckering  main  shock  and  is  presented  in  Sec¬ 
tion  5. 

Section  2  describes  another  study  of  near-regional  wave  forms 
from  small  foreshocks  and  aftershocks  associated  with  the  1967 
Koyna ,  India,  earthquake.  In  this  case,  a  different  distance 
range  as  well  as  crustal  structure  was  considered  with  the  same 
goals  in  mind.  Synthetic  seismograms  were  computed  to  gain  an 
understanding  of  the  wave  propagation  at  distances  of  100  to  150 
km  in  order  to  investigate  possible  depth-dependent  phases  in  the 
observed  data.  This  study  was  less  successful  because  crustal 
triplications  generally  obscured  the  interference  between  direct 
and  surface  reflections  in  the  body  waves.  General  inferences  on 
depth  were  possible  if  short-period  Rg  was  observed  in  the  data. 

The  general  conclusion  of  these  studies  on  short-period  wave 
propagation  was  that  short-period  wave  forms  can  be  interpreted 
to  infer  source  depth  accurately  if  an  appropriate  distance  range 
is  available  for  the  data  and  if  wave  propagation  effects  due  to 
the  crustal  model  are  understood.  This  research  will  be  extended 
in  future  work  to  consider  data  collected  in  other  areas  and  at 
other  distance  ranges.  It  is  conceivable  that  isolated  phases, 
such  as  Pn ,  may  be  usefully  modeled  at  short-periods  to  yield 


accurate  depth  estimates. 

Sections  3  and  4  are  concerned  with  the  theoretical  develop¬ 
ment  of  a  method  to  compute  synthetic  seismograms  in  layered 
media  with  corrugated  boundaries.  A  layer-over-half space  model 
was  considered  with  the  free  surface  of  the  layer  allowed  to  be 
corrugated.  SH  and  P-SV  line  sources  were  considered  in  Sections 
3  and  4,  respectively.  This  kind  of  problem  is  important  in  con¬ 
sidering  the  effect  of  waveguide  imperfections  on  the  propagation 
of  regional  seismic  phases  of  all  types.  The  results  show  that 
even  simple  periodic  corrugations  cause  very  complex  spectra  for 
elastic  waves.  Displacements  were  considered  on  the  free  surface 
as  well  as  in  the  lower  half space. 

A  reconnaissance  study  was  made  of  earthquake  sources  in  the 
African  continent  and  is  presented  in  Section  6.  The  purpose  of 
this  study  was  to  determine  accurate  earthquake  source  mechanisms 
and  depths  to  use  in  analysis  of  regional  wave  forms  from  the 
same  events.  Although  CTBT  issues  are  mainly  concerned  with 
regional  wave  propagation  in  the  Soviet  Union,  it  is  important  to 
gain  experience  in  effects  of  wave  propagation  in  different  crus¬ 
tal  and  mantle  environments.  Shield  areas,  for  example,  have 
certain  similarities  in  Q  and  velocity  structures.  These  simi¬ 
larities  are  stressed  when  it  is  important  to  infer  the  character 
of  wave  propagation  in  areas  of  the  Soviet  Union  for  which  there 
is  little  data  from  the  character  of  wave  propagation  in  other, 
known  areas.  The  African  continent  is  an  area  which  has  not  been 
studied  in  detail  but  has  a  nuclear  test  site  as  well  as  a  net¬ 
work  of  stations  and  regional  earthquakes.  The  results  of  Sec¬ 
tion  6  are  being  used  in  an  on-going  study  of  the  depth  effect  in 
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regional  seismograms. 

Section  7  contains  a  study  of  upper  mantle  structure  as 
determined  by  P  wave  forms.  This  work  was  supported  by  a  previ¬ 
ous  AFOSR  contract  and  carried  over  into  the  present  AFGL  con¬ 
tract.  A  basic  premise  in  studies  of  body  wave  propagation  in 
the  upper  mantle  is  the  effect  of  the  seismic  source  is  suffi¬ 
ciently  well  known  to  infer  the  smaller  effects  of  triplications 
due  to  mantle  discontinuities.  In  this  study,  an  intermediate 
depth  earthquake  which  occurred  beneath  Hispaniola  was  used  to 
examine  the  effect  of  mantle  triplications  for  Eastern  U.S.  and 
Western  Atlantic  structure.  It  was  thought  that  the  effective 
source  function  of  the  intermediate  depth  event  would  be  simpler 
than  those  used  in  previous  studies  since  near-source  reflections 
and  reverberations  could  be  avoided.  In  this  way  the  mantle  tri¬ 
plication  effect  would  be  more  clearly  seen.  The  results  show 
that  mantle  effects  are  still  difficult  to  model  in  the  long- 
period  P  wave  data  even  for  this  event.  However,  intermediate- 
period  data  recorded  at  DWWSSN  stations  clearly  showed  the  tri¬ 
plication  from  the  670  km  discontinuity.  Such  mantle  studies 
would  benefit  from  data  recorded  between  standard  long-period  and 
short-period  frequency  bands. 
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Abstract 


Wave  forms  for  eleven  foreshocks  and  48  aftershocks  of  the  Ms 
6.8  Meckering  earthquake  recorded  at  the  WWSSN  station  MUN  are 
analyzed  to  determine  the  depth  distribution  of  faulting  during 
this  unusual  intraplate  earthquake  sequence.  Clear  depth  phases 
including  sP  and  Rg  are  seen  in  the  local  seismograms  at  dis¬ 
tances  of  60  to  95  km  and  are  studied  using  synthetic  seismograms 
computed  using  generalized  ray  theory  and  wavenumber  integration 
techniques.  Large  Rg/S  ratios  seen  on  the  vertical  component 
short-period  data  for  many  events  imply  source  depths  less  than  2 
km.  The  short-period  P  wave  form  contains  the  best  depth  estima¬ 
tor  in  the  form  of  sP  so  that  depth  can  be  estimated  to  within  an 
uncertainty  of  about  1  km  for  most  events.  The  foreshocks  clus¬ 
ter  at  less  than  2  km  depth  and  most  aftershocks  occur  within  1 
km  of  the  surface.  A  few  aftershocks  occur  as  deep  as  7  km. 

These  results  are  consistent  with  a  previous  teleseismic  wave 
form  study  in  which  faulting  was  inferred  to  start  near  the  sur¬ 
face  at  1 . 5  km  depth  with  rupture  proceeding  downward  and  not 
exceeding  6  km  depth.  These  results  coupled  with  previous  stress 
studies  in  the  Australian  shield  and  models  of  crustal  strength 
show  that  faulting  in  a  cold  shield  area  is  a  near-surface  phe¬ 
nomenon  and  implies  that  most  of  the  crust  is  too  strong  to  be 
fractured . 


V> yy >  v  V  V  ■>  JV  *Y /  .  /  >  ■  v  r ,  - .  x  *  f-,, 


■vw 

V-S. 


9 


tralian  shield  (Everingham  and  Doyle,  1969,  Gordon,  1971;  Gordon 
and  Lewis,  1980).  The  main  event  was  the  subject  of  r.  teleseis- 
mic  body  wave  study  by  Vogfjord  and  Langston  (1987).  Also  incor¬ 
porating  a  static  dislocation  model  of  published  releveling 
measurements  across  the  fault  scarp,  they  determined  that  fault¬ 
ing  during  the  main  event  did  not  exceed  a  depth  of  about  six  km. 
Interpi etation  short-period  P  uaves  also  suggested  tha^  fault- 
ing  initiated  near  the  surface  at  about  1.5  km  depth  and  propa¬ 
gated  downward  to  a  similar  depth. 

The  unusual  nature  of  extremely  shallow  rupture  initiation, 
surface  faulting,  and  limited  depth  extent  of  faulting  was 
explained  by  a  strength  model  of  the  crust  after  Meisner  and 
Strehlau  (1982)  and  Strehlau  (1986).  Measured  heat  flow  in  the 
area  is  about  40  mW/m2  (Cull  and  Denham,  1979)  suggesting  that 
the  crust  is  relatively  cold  and  should  fail  in  a  brittle  fashion 
following  a  frictional  sliding  strength  relationship  down  to  Moho 
depths.  The  brittle-ductile  transition  is  substantially 
depressed  compared  to  a  crustal  column  with  higher  heat  flow, 
which  effectively  makes  the  crustal  strength  almost  everywhere 
much  higher  than  plausible  ambient  stresses.  Near  surface  stress 
measurements  in  the  area  (Denham  et  al,  1980)  show  relatively 
high  deviatoric  stress  of  about  20  MPa  near  the  surface  with  an 
assumed  depth  gradient  less  than  the  strength  curve  implied  by 
frictional  sliding.  The  only  area  within  the  crust  that  can 
fail,  according  to  this  model,  is  within  a  few  km  of  the  surface 
where  deviatoric  stress  exceeds  the  crustal  strength.  This  model 
is  significantly  different  from  models  proposed  for  interplate 
faulting  where  aseismic  slip  on  a  fault  in  the  brittle-ductile 
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dated  with  the  Meckering  event  and  wave  paths  are  entirely 
within  the  Archean  Yilgarn  block  of  Australia.  Immediately  to 
the  west  of  Mundaring  is  the  Perth  basin  bounded  by  the  Darling 
fault  which  underwent  several  kilometers  of  normal  offset  from 
Permian  time  throughout  the  Mesozoic.  See  Gordon  and  Lewis (1980) 
for  an  overview  of  geologic  structure  and  history  of  the?  area. 

The  structure  of  this  part  of  the  Yilgarn  block  has  been  stu¬ 
died  by  seismic  reflection  and  refraction  techniques  and  by  gra¬ 
vity  interpretation  (Mathur,  1974;  Drummond  and  Mohamed, 1986 ) . 
Table  1  displays  crustal  velocity  structures  determined  from 
these  studies  which  will  be  incorporated  in  synthetic  seismogram 
calculations  later  in  the  paper.  The  crust  between  the  Meckering 
area  and  Mundaring  appears  to  be  relatively  homogeneous  with  no 
significant  sedimentary  section  near  the  surface.  Both  struc¬ 
tures  show  relatively  high  P-velocities  of  about  6.1  km/sec  near 
the  surface  indicative  of  the  crystalline  granitic  terrain  seen 
in  the  surface  geology.  There  is  a  substantial  intermediate 
velocity  layer  at  10  km  depth  in  the  Mathur (1974)  model  which 
grades  into  a  high  velocity  lower  crust  at  only  17  km  depth.  The 
more  recent  Drummond  and  Mohamed( 1 986 )  model  is  smoother  showing 
very  small  positive  velocity  gradients  with  depth  and  with  an 
area  of  higher  gradient  near  13  km.  Topography  between  Mundaring 
and  Meckering  shows  very  little  relief. 

A  list  of  foreshocks  and  aftershocks  above  Richter  magnitude 
3.0  for  the  Meckering  event  was  published  in  Gordon  and  Lewis 
(1980)  and  was  used  to  initially  examine  the  WV7SSN  archives  for 
data  recorded  at  MUN .  An  effort  was  made  to  obtain  all  principal 
foreshocks  and  major  aftershocks  which  occurred  within  three  days 


after  the  main  event.  Other  events  were  obtained  as  available. 
In  addition  to  events  listed  in  Gordon  and  Lewis(1980),  a  sig¬ 
nificant  number  of  somewhat  smaller  events  were  discovered  and 
used  in  the  analysis. 

Events  used  in  this  study  are  listed  in  Table  2  by  P  arrival 
time  at  MUN.  The  S-P  time  and  maximum  trace  amplitude  on  the 
vertical  component  were  read  where  possible  to  estimate  distance 
and  local  magnitude,  respectively.  Distance  was  computed  using 
phase  identifications  from  the  synthetic  seismogram  computations 
Local  magnitude  was  computed  using  the  Richter ( 1958 )  formula  for 
California  events  and  corrected  for  the  WWSSN  short-period 
instrument  magnification  at  MUN.  Magnitudes  computed  this  way 
generally  agreed  with  those  listed  in  Gordon  and  Lewis (1980)  to 
within  0.1  magnitude  units,  although  Doyle  et  al(1968)  state 
that  local  magnitude  at  MUN  is  computed  using  the  maximum  hori¬ 
zontal  motions.  The  size  of  events  examined  in  this  study  range 
from  about  Ml  2.5  to  4.2. 

Locations  for  those  events  shown  in  Figure  1  were  determined 
from  data  from  a  sparse  regional  network.  Location  error  is 
probably  better  than  25  km  (Doyle  et  al,  1968).  Distance  deter¬ 
minations  used  in  this  study  relied  on  the  S-P  time  and  compari¬ 
son  with  synthetics  rather  than  the  published  locations.  All 
observed  events  at  MUN  could  then  be  treated  in  a  consistent  man¬ 
ner.  S-P  times  could  be  read  to  at  least  0.5  sec  in  all  cases 
which  implies  a  distance  error  of  +4.5  km  for  either  crustal 
model  assuming  a  Vs/Vp  ratio  appropriate  for  a  Poisson  solid. 

Figure  2  shows  three-component  wave  forms  for  two  aftershocks 
which  occurred  on  14  October.  These  data  were  hand  digitized 


from  copies  of  the  analog  seismogram  from  MUN  and  accurately 
show  the  main  features  of  the  original  seismograms.  The  two 
events  shown  are  characteristic  of  the  foreshock  and  aftershock 
data.  A  large  number  of  events  show  close  similarity  to  the  7-J8 
aftershock  in  which  the  vertical  wave  form  contains  a  relatively 
small  P  arrival  and  a  major  secondary  phase  at  about  12  seconds 
after  P.  A  careful  comparison  of  the  horizontal  components  with 
the  vertical  shows  that  the  S  arrival  is  always  sbout  9  seconds 
after  P  ~ince  most  events  are  at  similar  distances.  The  large 
secondary  arrival  on  the  vertical  component  characteristically 
shows  normal  dispersion  and  arrives  at  the  expected  Rayleigh  wave 
time  for  the  crustal  structure.  It  is  very  difficult  to  examine 
the  relative  phase  of  peaks  and  troughs  of  this  Rayleigh  wave 
between  the  vertical  and  horizontal  components  since  the  data  are 
band  limited  by  the  WWSSN  short-period  instrument  response.  The 
data  are  also  difficult  to  digitize  because  of  typically  high 
Rayleigh  amplitudes  for  most  events  above  local  magnitude  3.0. 
However,  a  careful  comparison  of  the  vertical  and  EW  components 
for  this  presumed  Rayleigh  phase  seen  in  relatively  small  events 
shows  retrograde  elliptical  motion.  This  Rayleigh  phase  is  desig 
nated  Rg  in  Figure  2  after  Ewing  et  al  (1957).  They  use  the  name 
Rg  for  the  f undcanental  Rayleigh  wave  with  periods  of  8  to  12  sec¬ 
onds  propagating  within  continents  at  regional  distances.  The  Rg 
phase  as  used  here  is  understood  to  be  the  fundamental  Rayleigh 
wave  at  one  second  period  seen  at  relatively  short  distance. 

In  comparison  to  the  very  simple  wave  form  of  the  7:38  event 
is  the  8:38  event  (Figure  2)  which  shows  larger  body  wave  phases 
relative  to  the  expected  Rg  phase  and  greater  overall  complexity. 


The  P  arrival  appears  to  have  at  least  two  distinct  arrivals 
within  it  compared  to  the  7:38  P  wave  and  it  is  difficult  or 
impossible  to  see  the  dispersion  of  Rg  in  the  S  wave  coda.  Events 
of  this  type  are  more  rare  than  those  displaying  the  large  Rg 
phase . 

Figure  3  displays  six  aftershocks  showing  the  typical  variab¬ 
ility  seen  in  the  vertical  wave  form.  Usually  the  Rg  phase  is 
well  developed  but  may  range  in  size  from  being  comparable  to  the 
S  wave  amplitude  to  an  order  of  magnitude  larger  than  the  S  wave. 
Presumably ,  these  relative  amplitude  variations  are  caused  by 
different  source  mechanisms  between  events  and  source  depth.  The 
radiation  pattern  of  each  phase  will  depend  on  mechanism  and  the 
excitation  of  the  Rg  wave  should  be  expected  to  decrease  exponen¬ 
tially  with  source  depth.  The  initial  working  hypothesis  that  is 
suggested  by  this  cursory  analysis  of  the  data  is  that  it  may  be 
possible  to  estimate  bounds  on  the  source  depth  from  Rg  relative 
amplitudes  alone.  The  mitigating  factors  will  be  wave  propagation 
effects  for  the  body  wave  phases  in  the  crust  and  the  unknown 
source  mechanism. 

Only  the  vertical  wave  forms  were  used  in  the  analysis  that 
follows.  Since  these  events  are  not  accurately  located  and 
because  of  the  practical  difficulties  of  digitizing  short-period 
data  from  analog  records,  the  horizontal  components  of 

motion  was  confined  to  determ  arrival  times.  Any  vector 

rotation  of  these  horizontal  data  to  obtain  radial  and  tangential 
motions  for  study  would  probably  be  of  very  poor  quality  because 
of  unknowns  in  back  azimuth  to  the  source  due  to  epicenter  error 
and  because  of  the  digitization.  In  fact,  many  of  the  Rg  phases 


shown  in  the  figures  were  obtained  by  digitizing  the  peaks  and 


troughs,  which  could  be  seen  on  the  records,  and  then  linearly 


interpolating  for  plotting  purposes.  The  Rg  analysis  of  all 


events  was  further  simplified  by  simply  measuring  the  peak  to 


peak  amplitudes  directly  from  the  seismogram.  Likewise,  P  and  S 


amplitudes  were  also  read  although  the  P  wave  form  is  used  for 


another  purpose  below. 


Synthetic  Seismogram  Calculation 


Synthetic  seismograms  were  computed  to  investigate  the  exci 


tation  of  the  Rg  phase  relative  to  the  body  wave  phases  as  a 


function  of  source  depth.  Because  the  total  wave  field  is  needed 


for  this  problem,  a  wavenumber  integration  method  was  employed  to 


compute  the  Green’s  functions  required  for  constructing  moment 
tensor  point  sources.  Details  of  the  technique  may  be  found  in 
Barker ( 1984 )  and  Apsel(1979).  An  integral  over  real  wavenumber 
is  performed  using  an  adaptive  Filon  method.  The  kernel  of  the 
integral  contains  a  product  of  the  vertical  wave  function  and  a 


Bessel  function  of  integer  order  depending  on  the  type  of  seismic 


source.  The  vertical  wave  function  is  computed  using  a  propaga¬ 


tor  matrix  formalism  for  plane  elastic  layers  which  are  homoge¬ 


neous  and  isotropic.  Anelastic  attenuation  is  included  using 


complex-valued  wave  velocities.  The  moment  tensor  formalism 
given  by  Langston( 1981 ,  1983)  is  used  to  combine  the  four 


Green’s  functions  needed  for  computation  of  vertical  displacement 


from  a  general  moment  tensor  source. 


Vertical  displacement  seismograms  for  ranges  of  60  to  90  km 


hhsm 


were  calculated  for  the  Mathur(1974)  model  of  Table  1.  The 
short-period  WWSSN  instrument  response  and  a  Brune  (1970)  time 
function  with  a  corner  frequency  of  1  hz  were  convolved  with 
these  displacement  Green’s  functions  to  produce  the  final  syn¬ 
thetic  seismograms.  The  wavenumber  integration  was  carried  out 
over  the  interval  between  0.01  and  0.4  sec/km  (in  ray  parameter) 
which  was  sufficient  to  include  all  body  and  surface  wave  phases 
for  the  structure  model.  A  Nyquist  frequency  of  5  hz  was  assumed 
and  256  frequency  samples  computed  to  yield  a  time  series  after 
inverse  Fourier  transforming  of  512  points  at  0.1  sec  sampling. 

Figure  4  displays  vertical  displacement  synthetic  seismograms 
for  a  range  of  80  km,  the  average  distance  of  the  foreshocks  and 
aftershocks  from  MUN.  Wave  forms  for  five  different  source 
depths  show  that  the  excitation  of  the  Rg  wave  decreases  with 
depth  as  expected,  relative  to  the  body  wave  phases,  and  is  only 
well  developed  for  sources  above  4  km  depth.  The  character  of 
the  wave  forms  for  the  0.25  km  source  depth  is  similar  to  some  of 
the  data  wave  forms  of  Figure  3. 

At  this  point  it  should  be  mentioned  that  the  original  Mathur 
(1974)  crustal  model  did  not  have  the  thin  low  velocity  layer  at 
the  top  of  the  section  as  shown  in  Table  1.  Initial  synthetic 
seismogram  calculations  without  the  layer  produced  a  single 
undispersed  Rayleigh  pulse.  The  published  refraction  data  did 
not  resolve  the  upper  0.5  km  of  the  crust.  However,  the  disper¬ 
sion  of  the  short-period  Rg  data  clearly  suggests  that  a  posi¬ 
tive  velocity  gradient  exists  near  the  surface.  Theoretical  dis¬ 
persion  curves  given  by  Ewing  et  al(1957)  were  used  to  approxi¬ 
mately  match  the  observed  dispersion  and  to  obtain  an  equivalent 


layer  thickness  for  a  layer  with  a  modest  velocity  contrast 
appropriate  for  near  surface  granitic  rocks  (Press,  1966).  It 
was  important  to  produce  a  dispersed  Rg  wave  form  more  in  accord 
with  the  data  to  investigate  the  amplitude  dependence  of  the  Rg 
wave  with  near  surface  structure.  Dispersion  will  tend  to  lower 
maximum  amplitudes  since  energy  is  spread  out  over  time.  How¬ 
ever,  the  net  effect  of  the  lower  velocities  near  the  surface  was 
to  increase  amplitude  in  nearly  the  same  amount  that  dispersion 
decreased  amplitude,  in  comparison  to  the  body  waves.  The  behav¬ 
ior  of  the  body  waves  remained  the  same  between  models. 

Figure  5  shows  the  typical  variation  of  relative  amplitude 
expected  as  a  functioi  of  range.  There  are  some  changes  in 
relative  amplitudes  mainly  due  to  structure  effects  on  the  body 
waves.  In  fact,  maximum  amplitudes  of  the  wave  form  appear 
roughly  constant  over  this  distance  range  because  of  wide  angle 
reflections  from  the  shallow  intermediate  crustal  layer.  This 
overcomes  the  differences  in  geometric  spreading  expected  by  a 
simple  direct  body  wave  (1/r)  and  a  simple  surface  wave  (l/r°-5). 

It  can  be  expected  that  variations  in  Q  and  scattering  along 
the  path  of  transmission  will  affect  Rg  amplitudes.  The  Q  struc¬ 
ture  assumed  in  the  seismogram  calculation  is  reasonable  for  a 
continental  area  but  may  be  somewhat  high.  However,  a  1  hz  sur¬ 
face  wave  propagating  over  80  km  is  only  attenuated  to  about  80% 
of  its  nominal  value  with  this  structure.  Using  a  Q  of  250 
instead  of  500  for  shear  waves  reduces  a  1  hz  surface  wave  to  70% 
of  its  unattenuated  value.  Scattering  is  probably  minimal  along 
the  path  since  topographic  changes  are  minimal.  Most  effects  of 
scattering,  if  they  occur,  will  tend  to  reduce  surface  wave 


amplitudes  while  lengthening  the  surface  wave  coda  with  scattered 
energy.  Thus,  the  seismograms  presented  here  will  maximize  the 
expected  Rg  phase  amplitudes.  Any  effects  of  Q  and  scattering 
will  reduce  the  observed  Rg  amplitudes  tending  to  make  an  event 
appear  deeper  than  it  really  may  be.  Therefore,  source  depths 
interpreted  from  these  synthetic  Rg  amplitudes  will  tend  to  be 
overestimated,  all  other  factors  being  equal. 

A  straightforward  interpretation  of  Figure  4  would  suggest 
that  Rg  relative  amplitudes  of  one  or  more  are  only  possible  for 
sources  less  than  4  km  in  depth.  However,  it  is  possible  that 
source  orientation  may  affect  the  amplitudes  of  all  phases  such 
that  Rg  may  be  near  a  radiation  pattern  minimum  when  the  body 
waves  are  large,  even  for  a  very  shallow  source.  To  investigate 
this,  a  grid  search  process  was  set  up  to  sample  all  possible 
fault  orientations  for  Rg/S  and  Rg/P  ratios.  A  10  degree  incre¬ 
ment  was  assumed  for  changes  in  dislocation  dip,  rake,  and 
strike.  The  maximum  Rg,  S,  and  P  amplitudes  were  found  for  spe¬ 
cified  time  windows  of  the  synthetic  seismogram  which  included 
these  phases,  ratios  taken,  and  the  ratios  were  then  placed  in 
bins  with  an  integer  interval.  It  was  found  that  the  Rg/S  ratio 
was  the  most  robust  measure  of  source  depth  since  it  appeared 
possible  to  get  nearly  any  Rg/P  ratio  for  any  source  depth.  This 
is  logical  since  the  Rayleigh  wave  from  a  dislocation  source  is 
closely  tied  to  the  shear  wave  potential  and  it  is  not  difficult 
to  produce  nodal  P  waves  for  specified  azimuths  and  incidence 
angles . 

For  a  source  depth  of  0.25  km  the  range  of  Rg/S  varied  from  4 
to  20.  The  range  of  values  for  Rg/S  assuming  other  source  depths 


varied  generally  from  0  to  a  well  defined  maximum.  Figure  6 
shows  the  maximum  Rg/S  ratio  as  a  function  of  source  depth  and 
range  for  these  trials.  This  graph  is  to  be  interpreted  as  pro¬ 
viding  a  maximum  estimate  for  source  depth  based  on  a  measured 
Rg/S  ratio.  High  Rg/S  ratios  of  10  or  greater  imply  sources  with 
depths  less  than  2  km  for  observations  within  the  assume  distance 
ranges.  A  low  Rg/S  ratio  may  still  be  possible  for  a  2  km  depth 
source  but  there  is  no  way  from  the  ratio  itself  to  determine 
this.  The  low  ratio  will  be  interpreted  as  being  due  to  a  larger 
source  depth  with  the  understanding  that  this  will  be  maximum 
estimate  of  depth. 

In  addition  to  the  Rg  wave  depth  effect,  note  the  character 
of  the  body  wave  phases  in  Figure  4.  As  the  source  becomes 
deeper,  the  simple  P  pulse  for  near-surface  sources  splits  into  a 
distinctive  double  arrival.  The  isotropic  source  does  not  show 
this  effect,  suggesting  that  the  second  arrival  in  the  other 
synthetics  is  an  S  to  P  conversion.  Indeed,  an  analysis  of  the  P 
wave  forms  using  ray  and  generalized  ray  theory  shows  that  the 
secondary  phase  after  direct  P  is  the  surface  reflection  sP. 

Body  wave  phases  may  be  expected  to  be  sensitive  to  details 
of  structure.  In  the  Mathur(1974)  model,  the  initial  P  arrivals 
are  dominated  by  wide  angle  reflections  from  both  top  and  bottom 
interfaces  of  the  6.7  km/sec  layer  (Table  1).  Alternatively,  the 
Drummond  and  Mohamed( 1986)  model  consists  of  smooth  gradient 
zones  throughout  the  crust  and  would  produce  true  turning  rays. 
Figure  7  displays  a  synthetic  P  wave  record  section  for  a  CLVD 
source  at  4  km  depth  incorporating  the  velocity  structure  of  both 
crustal  models.  These  synthetics  were  computed  used  generalized 


ray  theory  (Helmberger  and  Harkrider,  1978).  The  gradient  model 
was  constructed  by  using  a  number  of  thin  layers  and  responses 

computed  by  summing  all  primary  reflections  from  each  interface 
in  the  model. 

Since  the  Mathur  model  consists  of  only  a  few  discrete 
layers,  the  wave  forms  show  several  small  arrivals  between  P  and 
sP  due  to  wide  angle  pP  reflections  from  the  two  interfaces  of 
the  6.7  km/sec  layer  (column  A,  Figure  7).  The  gradient  model 
results  (column  B)  display  simpler  wave  forms  overall  but  contain 
the  emergence  of  a  phase  at  distances  greater  than  90  km.  The 
high  velocity  gradient  near  the  Moho  produces  a  triplication  in 
the  wave  form  for  these  distances.  Even  with  this  complication, 
sP  is  still  an  easily  seen  phase  in  the  wave  form.  Overall,  both 
crustal  models  allow  for  the  development  of  an  unambiguous  sP 
phase.  However,  it  is  clear  that  relative  amplitudes  will  depend 
on  model  details.  The  phase  pP  is  not  well  developed  at  this 
range  primarily  because  of  the  inefficiency  of  the  free  surface 
reflection  at  the  appropriate  ray  parameters.  The  S-to-P  conver¬ 
sion  is  near  a  maximum  while  P  reflection  is  near  a  minimum. 

Theoe  modeling  results  suggest  that  secondary  phases  i„i  the  P 
wave  form,  if  seen,  will  be  sP.  Source  depth  can  be  determined 
from  such  arrivals  by  simply  reading  the  sP-P  time  and  using  the 
ray  theory  approximation 


lsP-P  ~H{%  +  r)p) 


to  solve  for  depth.  Here  0Q  and  ^  are  the  vertical  slow- 


nesses  for  P  and  S  waves, 


respectively.  The  appropriate  ray.  par 
aneter,  p,  is  1/6.2  sec/km.  This  equation  was  tested  by  using 
the  wavenumber  synthetic  seismograms  (e.g.  Figure  4).  The  sP-P 
times  were  read  from  the  synthetic  wave  forms,  depth  calculated 
from  equation  (1)  and  these  values  compared  to  the  actual  depths 
used  in  the  synthetic  calculations.  The  depth  determined  by  ray 

theory  agreed  to  within  0.5  km. 

If  sP  phases  can  be  seen,  then  it  is  possible  to  infer  abso¬ 
lute  source  depth  very  accurately.  Timing  between  peaks  on 
observed  wave  forms  can  easily  be  read  to  0.3  sec  which  gives  a 
possible  depth  error  of  only  1  km.  This  will  be  taken  to  be  the 
error  of  the  determination  if  the  separation  between  phases  is 
clear.  However,  there  still  exists  the  possibility  that  other 
crustal  structure  arrivals,  like  the  arrivals  from  the  Moho  of 
the  Drummond  and  Mohamed  model  shown  in  Figure  7,  may  make  this 
ambiguous.  In  this  case,  we  rely  on  the  data  to  make  an  argument 
in  self-consistency.  It  was  seen  that  events  with  the  same  S-P 
times  can  have  a  single  P  arrival,  implying  a  very  shallow  source 
depth,  or  two  clear  P  arrivals  implying  P  and  sP  and  a  larger 
depth.  Also,  it  was  seen  that  events  with  a  single  P  arrival 
occurred  at  ranges  from  68  to  89  km  implying  that  the  Drummond 
and  Mohamed  model  may  need  some  modification  in  structure  near 
the  Moho  to  push  the  triplication  out  to  greater  ranges. 


The  experience  gained  by  examining  the  synthetic  seismograms 
suggests  two  ways  of  treating  the  data  to  determine  source  depth. 
The  first  consists  of  estimating  the  maximum  allowable  source 
depth  from  the  Rg/S  ratio.  To  do  this,  it  is  only  necessary  to 
read  the  maximum  amplitudes  of  the  Rg  and  S  phases  from  seismo¬ 
gram  copies,  take  the  ratio,  and  read  off  the  maximum  depth  from 
the  curves  of  figure  6.  The  other  method  consists  of  examining 
the  P  arrivals  to  determine  the  timing  of  inferred  sP  phases  and 
compute  the  source  depth  from  equation  (1).  These  ax*e  very 
simple  techniques  for  interpreting  the  data  but  are  probably  the 
only  techniques  warranted  by  the  data  at  this  time.  The  use  of 
wave  form  inversion  techniques  to  determine  mechanism  or  other 
parameters  depends  the  veracity  of  the  crustal  model  used  and 
the  quality  of  the  digital  wave  form  data.  Since  the  data  are 
from  a  single  station,  are  not  conducive  to  accurate  digitiza¬ 
tion,  and  because  crustal  models  for  the  area  are  still  only 
approximate  it  seems  unwise  to  pursue  such  inversions. 

Vert.vcal  component  seismograms  for  the  events  listed  in  Table 
2  were  examined  to  obtain  the  Rg/S  ratio.  This  ratio  was  taken 
for  those  events  where  the  Rg  phase  was  on  scale  and  is  listed  in 
the  table.  P  waves  of  events  with  clear  P  phases  were  digitized 
for  display  purposes  to  show  the  picks  of  sP.  The  sP-P  time  was 
read  from  the  digitized  wave  form  and  interpreted  to  obtain 
source  depth.  All  depth  results  are  shown  in  Table  2. 

Figure  8  shows  vertical  wave  forms  of  foreshocks  which  were 
clear  enough  on  the  record  to  be  digitized.  Note  that  the  wave 


A  compilation  of  maximum  source  depths  (Figure  14)  inferred 
from  Rg/S  ratios  shows  about  17  events  with  defau.it  depi  hs  at 
to  6  km  but  with  most  events  occurring  above  3  km.  Only  the  8:04 
event  of  Figure  11  showed  no  obvious  Rg  wave  which  could  be 
interpreted.  This  aftershock  was  inferred  to  be  at  7.6  km  depth 
from  it’s  sP-P  time  and  is  the  deepest  event  found.  The 
frequency  of  events  with  source  depth  found  using  sP  times  is 
shown  in  Figure  15.  Again,  most  of  the  events  occurred  above  3 
km  but  it  appears  that  seismicity  is  concentrated  near  the  sur¬ 
face  and  falls  off  exponentially  with  depth. 

Figure  14  is  biased  towards  ‘deep’  hypocenters  because  of  the 
nature  of  the  Rg  analysis.  Figure  15  may  be  slightly  biased 
towards  very  shallow  events  for  those  interpretations  which  may 
be  affected  by  nodal  sP.  There  may  also  be  some  bias  in  the 
choice  of  data  in  this  magnitude  range  but  it  is  not  known  in 
what  direction  the  bias  could  push  the  distribution.  A  study  of 
the  records  showed  many  smaller  events  which  could  be  seen  only 
by  their  Rg  waves.  If  these  events  were  included  they  would  cer¬ 
tainly  increase  the  number  of  events  at  depths  shallower  than  3 
km.  However,  events  of  the  same  moment  that  occurred  deeper 
would  be  undetectable  since  they  would  not  generate  high  ampli¬ 
tude  Rg  phases.  In  the  largest  events  Rg  went  off  scale  and  was 
unavailable;  these  are  the  events  in  Table  2  for  which  Ml  is  not 
given.  However,  these  larger  events  concentrate  at  about  2  km 
depth  as  determined  by  sP.  Thus,  it  appears  that  very  shallow 
source  depths  are  characteristic  of  this  earthquake  sequence. 


Gurrentiy  unique  among  'the  seismogemc  regi 
Some  models  of  faulting,  particularly  a 
areas  such  as  the  San  Andreas  fault  zone,  h 
affected  by  ductile  rheology  below  the  brit 
region.  Slip  in  the  ductile  region  before 
the  crack  tip  near  the  brittle-ductile  tran 
instability  for  the  occurrance  of  future  ev 
(Thatcher,  1983;  Strehlau,  1986).  An  alter 
occurring  during  an  earthquake  in  the  upper 
large  stresses  to  build  up  at  the  brittle-d 
causing  the  same  kind  of  instability.  Seisi 
almost  exclusively  in  the  upper  brittle  lay* 
rupture  into  the  ductile  region  (Strehlau,  ; 
ates  in  this  high  strength  transition  zone  1 
and  ductile  layers  of  the  lithosphere.  A  v< 
model  has  been  proposed  by  Gordon  and  Wellme 
and  Lewis (1980)  for  the  occurrence  of  the  Me 
suggest  that  aseismic  slip  on  a  deep  shear  z 
stressed  the  upper  portion  of  the  crust  to  e 
ure  during  the  Meckering  earthquake. 

In  contrast  with  this  model  is  the  infer 
of  the  Meckering  event  which  initiated  near 


upper 


third  of  the  fault  plane  and  propagated  downward  to  termi¬ 
nate  at  only  6  km  depth.  This  behavior  is  consistent  with  a  model 
suggested  by  Furlong  (personal  communication)  and  Vogfjord  and 
Langston (1 987)  which  was  described  above.  In  this  model  the 
strength  of  the  Australian  shield  exceeds  the  ambient  deviatonc 
stresses  everywhere  except  near  the  surface.  The  fault  is  driven 
simply  by  high  deviatoric  stresses.  The  high  strength  of  the 
crust  requires  that  faulting  stops  before  penetrating  too  deeply 
into  the  crust.  Rupture  proceeds  down-dip  in  a  manner  suggested 
by  Strehlau (1986)  until  dynamic  frictional  resistance  exceeds 
the  stresses  induced  by  faulting  at  the  crack  tip. 

The  southwest  seismic  zone  of  Australia  may  be  an  area  that 
has  interesting  possibilities  for  physical  experiments.  Because 
the  hypocentral  region  is  only  at  2  km  depth  or  so  it  would  be 
very  feasible  to  drill  into  it  with  existing  technology.  Gauges 
could  be  emplaced  over  the  length  of  the  borehole  to  measure 
physical  properties  such  as  temperature,  strain,  stress,  and  the 
state  of  pore  fluids.  Such  an  area  would  be  an  accessible  earth¬ 
quake  laboratory  in  which  physical  and  seismic  monitoring  would 
yield  information  on  the  formation  of  continental  intraplate 
events.  Because  of  the  homogeneity  of  crustal  structure  and 
topography  in  the  area,  complications  of  structure  normally  found 
in  most  tectonic  areas  of  the  world  would  be  avoided  or  mini 
mized.  It  may  also  be  possible  to  perform  limited  physical 
experiments  in  the  generation  of  man-made  earthquakes  depending 
on  physical  conditions  found  in  the  upper  crust  (e.g.,  Rayleigh 
et  al,  1972).  In  any  case,  scientific  study  of  such  an  area 
would  be  complimentary  to  those  sites  along  active  interplate 


usefully  employed.  Indeed,  it  is  possible,  in  principle,  to 
employ  wave  form  analysis  in  those  areas  which  are  more  compli¬ 
cated  but  have  been  extensively  studied  so  that  wave  propagation 
effects  are  well  understood.  These  kinds  of  methods  would  sup¬ 
plement  network  analysis  of  local  seismicity  in  determining 
source  parameters  even  for  those  networks  which  are  instrumented 
with  only  short-period  vertical  seismometers. 
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Figure  1 .  Sketch  map  of  the  Meckering  area  showing  epicenters  of 
earthquakes  in  the  SW  Australian  seismic  zone,  the 
Meckering  fault  scarp,  and  location  of  Mundaring  (MUN) 
station.  Map  was  modified  after  Figure  1  of  Denham  et 
al  (1980). 

Figure  2:  Two  aftershocks  of  14  October  (see  Table  2)  recorded 
on  the  WWSSN  short-period  system  at  MUN.  Maximum 
ground  motion  for  the  vertical  (Z),  East-West  (E) ,  and 
North-South  (N)  components  is  normalized  and  given  in 
microns  for  each  trace.  Major  phases  are  denoted  by 
P,  S,  and  Rg.  Most  events  have  large  Rg  phases  evi¬ 
dent  on  the  vertical  components  as  shown  by  the  7 : 38 
event.  In  contrast,  the  8:38  event  has  a  poorly 
developed  Rg  phase  in  addition  to  a  complex  P  wave. 


Figure  3:  Short-period  vertical  component  wave  forms  from  six 
aftershocks  of  14  October  (Table  2)  showing  the 
variability  of  relative  amplitudes  among  the  major 
arrivals . 


Figure  4:  Synthetic  seismograms  for  the  four  fundamental  terms 
needed  to  construct  a  general  moment  tensor  point 
source.  \SS,  VDS,  CLVD,  and  ISO  denote  vertical 
strike-slip,  vertical  dip-slip,  compensated  linear 
vector  dipole,  and  isotropic  sources,  respectively. 
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Distance  to  the  receiver  is  80  km  and  seismograms  are 
shown  as  a  function  of  source  depth,  H.  The  wavenum¬ 
ber  integration  technique  was  used  calculate  these 
synthetics  for  the  Mathur(1974)  model  of  Table  1. 

They  include  the  short-period  WWSSN  instrument 
response  and  a  Brune(1970)  time  function  with  a  corner 
frequency  of  1  hz .  Amplitude  is  in  cm  for  an  assumed 
source  moment  of  1020  dyn  cm.  Note  the  existence  of 
large  Rg  phases  for  sources  at  depths  of  2  km  and 
less  compared  to  P  and  S  phases.  Rg  is  not  well 
developed  for  deepei  sources.  Also  note  the  existence 
of  sP  which  separates  from  the  first  arriving  P  as 
source  depth  increases.  The  isotropic  source  does  not 
show  sP  since  this  source  does  not  generate  S  waves. 

Figure  5:  The  effect  of  range  on  wave  forms.  A  vertical  strike- 
slip  ( VSS )  source  at  4  km  depth  is  assumed.  Model 
parameters  are  the  same  as  those  used  in  Figure  4 . 

Figure  6:  Maximum  Rg/S  ratio  plotted  as  a  function  of  source 

depth.  Curves  are  shown  for  ranges  of  60,  70,  80,  and 
90  km.  Note  that  large  Rg/S  ratios  of  10  or  more 
imply  shallow  source  depths  of  less  than  2  km.  How¬ 
ever,  small  Rg/S  ratios  may  also  be  produced  for 
sources  less  than  2  km  depth  because  of  source  radia¬ 
tion  pattern  effects.  These  curves  are  used  to  esti¬ 
mate  maximum  allowable  source  depths  from  the  Rg  data. 
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The  effect  of  crustal  model  on  expected  P  wave  forms. 
Column  A  shows  a  profile  of  P  waves  computed  using 
generalized  ray  theory  and  the  Mathur(1974)  model. 
Range  is  given  in  kilometers.  Source  depth  is  6  km. 

A  CLVD  source  is  assumed  since  it  is  likely  that  many 
Meckering  events  are  caused  by  thrusting;  a  CLVD  seis¬ 
mogram  can  be  interpreted  as  being  excited  by  a  45 
degree  dipping  thrust  fault  seen  at  45  degrees  azi¬ 
muth.  Column  B  shows  a  profile  of  synthetic  seismo¬ 
grams  again  computed  using  generalized  ray  theory  but 
assuming  the  Drummond  and  Mohamed ( 1986 )  model.  Note 
that  the  relative  amplitudes  of  sP  and  P  depend  on 
range  and  the  type  of  crustal  model  in  detail  but  sP 
is  always  prominent.  At  ranges  of  90  and  100  km  the 
Drummond  and  Mohamed  model  produces  a  triplication  due 
to  the  high  velocity  gradient  at  the  Moho.  This 
causes  prominent  Moho  turning  phases  like  sPm.  There 
is  little  evidence  for  this  triplication  in  the  data 
although  there  are  few  events  seen  at  distances 
greater  than  90  km. 

Observed  short-period  vertical  component  wave  forms 
for  six  foreshocks.  Amplitude  is  in  microns  for  maxi¬ 
mum  displacement.  The  large  excursions  in  the  data 
for  the  23:02  and  23:13  events  after  the  Rg  phases  are 
from  the  effect  of  minute  marks  which  were  digitized. 
Aside  from  different  signal  to  noise  ratios  which 
principally  depend  on  the  size  of  each  event,  note  the 
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to  determining  the  geological  relationships  and  perhaps  the 
cause  of  this  reservoir  induced  seismicity.  The  solution  to 
this  problem  is  important  because  the  more  we  understand 
about  induced  seismicity  the  better  we  can  avoid  the 
destruction  caused  by  these  "man-made''  events  in  the  future. 

Due  to  the  sparse  station  coverage,  the  shallow  source 
depth,  and  the  small  magnitude  of  most  of  the  events  in  the 
Koyna  area,  accurate  determination  of  detailed  source 
parameters  has  been  very  difficult.  In  this  study,  thirteen 
events  that  occurred  in  December  1967  are  examined  to  see  if 
some  of  the  problems  in  determining  accurate  source 
parameters  that  have  been  encountered  in  the  past  can  be 
overcome  through  the  study  of  the  regional  short-period  wave 
propagation.  If  aspects  of  the  regional  wave  propagation 
are  understood,  the  analysis  of  the  source  parameters  of  the 
smaller  events  may  be  possible.  In  order  to  do  this,  it 
must  be  determined  what  aspects  of  the  seismogram  are 
controlled  by  source  and  what  aspects  are  controlled  by 
earth  structure.  To  answer  this,  the  nature  of  the 
seismogram  at  local  and  regional  distances  and  the  factors 
that  control  its  character  must  be  addressed. 

The  nature  of  the  short-period,  regionally-recorded 
seismogram  is  usually  complicated  because  the  recorded 
waveforms  are  shaped  by  a  complex  interaction  of  seismic 
source  with  the  crustal  structure.  The  near- source,  near¬ 
receiver,  and  most  of  all,  the  regional  structure  between 
the  source  and  receiver  all  affect  the  character  of  the 
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seismogram.  It  is  these  complicated  effects  caused  by  the 
regional  variation  in  the  crustal  structure  wnich  usually 
outweigh  the  source  effects  and  control  the  character  of  the 
regie nal  seismogram.  Due  to  these  variations,  the  regional 
seismogram  is  made  up  of  a  variety  of  diffraction  effects 
such  as  head  waves,  critical  and  post-critical  reflections, 
and  surface  waves.  Unlike  long-period  studies,  where 
knowledge  of  the  gross  earth  structure  is  sufficient, 
short-period  studies  require  finer  structural  detail  to 
place  meaningful  constraints  on  the  source  parameters. 

The  events  that  are  studied  were  selected  because  the 
relative  simplicity  of  the  short-period  waveforms  suggested 
that  there  is  a  simple  earth  structure  that  could  be  easily 
modeled,  and  hence  source  details  could  be  inferred  from  the 
signal  information.  Even  though  the  events  are  small, 
approximately  magnitude  3.0,  and  insignificant  in  terms  of 
strain  and  energy  release,  they  provide  the  same  information 
as  the  larger  events  in  the  Koyna  area  as  to  the  nature  of 
faulting,  rupture  depth,  and  location.  Also,  since  there 
are  so  many  small  events,  the  redundancy  of  the  wave 
propagation  information  might  be  used  to  advantage  in 
identifying  specific  major  regional  wave  phases. 

In  this  study,  wave  propagation  effects  from  many  small 
events  recorded  at  a  single  station  are  examined.  The 
effectiveness  of  using  a  single  station  to  obtain  useful 
seismic  information  for  some  data  sets  has  been  demonstrated 
by  Langston  (1979,1987)  and  this  study  was  motivated  along 
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similar  lines.  The  earthquakes  were  recorded  at  the  station 
POO,  located  120  km  north  of  the  Koyna  Dam  at  Poona,  India. 
At  this  range,  POO  is  able  to  record  small  events  which  are 
not  recorded  at  stations  at  teleseismic  distances.  POO,  a 
WWSSN  station,  has  had  a  very  sensitive  short  period 
seismograph  in  operation  since  the  early  1950’s  (Guha  et 
al.,1971)  and  is  just  5°  off  azimuth  from  due  north  of  the 
dam  site.  Since  the  location  of  POO  is  almost  directly 
north  of  these  events,  the  horizontal  components  there  are 
nearly  naturally  polarized  into  radial  (north-south)  and 
tangential  (east-west)  directions. 

Synthetic  seismograms  utilizing  previously  determined 
velocity  structures  were  computed  using  two  different 
methods,  a  generalized  ray  theory  and  a  wavenumber 
integration  technique,  to  model  the  wave  propagation  and 
event  source  characteristics  of  the  events  in  the  Koyna  Dam 
area.  It  will  be  shown  that  with  knowledge  of  regional  wave 
propagation  even  short-period  data  recorded  at  a  single 
station  can  be  used  to  constrain  some  source  parameters  of 
small  earthquakes.  The  identification  of  the  major  regional 
body  wave  phases  and  the  existence  of  surface  waves  are  used 
to  constrain  source  depth  to  within  a  few  kilometers,  as 
well  as  infer  source  mechanism  and  event  distance  for  all 
the  events  studied.  The  simple  velocity  structure  and  its 
effect  on  wave  propagation  at  ranges  of  110-150  km  makes  the 
determination  of  these  source  parameters  possible. 

Th**  method  used  in  this  study  was  also  applied  to  an 


explosion  source  along  with  the  earthquake  sources.  In 
doing  this  it  was  determined  that  possibilities  exist  for 
the  discrimination  between  natural  events  and  nuclear  tests. 
Through  the  identification  of  the  source  parameters  the 
likelihood  of  a  nuclear  source  could  be  eliminated  by 
constraining  rupture  depths  to  >10  km  or  by  locating  events 
to  be  in  areas  unsuited  for  underground  testing  (a.g. 

Dahlman  and  Israelson , 1977 ) .  In  addition  the  method  can  be 
used  in  regions  where  limited  seismic  station  coverage 
renders  other  seismological  methods  of  identification 
ineffective . 

II.  DATA  PROCESSING 

Microfiche  copies  of  all  the  available  WWSSN  short-  and 
long- period  records  for  the  month  of  December  1967  recorded 
at  the  Poona,  India  station,  POO,  were  obtained  from  the 
National  Geophysical  and  Solar  Terrestrial  Data  Center 
EDS/NOAA  in  Boulder,  Colorado.  Contained  within  this  data 
set  were  both  foreshocks  and  aftershocks  of  the  December  10 
raainshock .  These  events  were  too  small  to  be  recorded  on 
the  long-period  system  so  only  short-period  seismograms  were 
studied . 

For  an  event  to  be  selected  for  digitization  it  had  to 
be  on  scale  and  be  free  of  signal  interference  from  adjacent 
events.  These  problems  were  encountered  many  times  since 
there  were  numerous  events  distributed  over  a  large  range  in 


magnitude.  The  events  that  were  finally  picked  for 
digitization  were  digitized  from  peak  to  peak  at  irregular 
intervals,  The  digitizing  was  done  in  this  manner  because 
the  short-period  seismogram  varies  quickly  from  peak  to  peak 
and  may  in  fact  only  be  visible  at  the  peak  of  the 
waveforms.  Other  difficulties  in  digitizing  included 
picking  out  the  peaks  of  the  F  waves  due  to  the  fact  that 
the  background  noise  was  nearly  the  same  size  as  the  P  wave 
amplitude  in  most  cases.  Also,  due  to  the  compressed  time 
scale  the  closeness  of  S  wave  peaks  made  it  hard  to  tell 
which  peak  was  next  in  the  continuous  signal. 

Once  digitized,  a  magnification  correction  was  applied 
to  the  data  and  it  was  interpolated  to  a  uniform  . 1  second 
interval  between  points.  The  final  interpolated  results  for 
the  three  components  of  the  thirteen  events  chosen  for  the 
study  are  illustrated  in  Figures  2.1,  2.2,  and  2.3. 

HI.  CRUSTAL  STRUCTURE  OF  KOYNA  REGION 

Because  of  the  previous  aseismic  nature  and  seemingly 
simple  geology  of  the  area,  the  geology  and  tectonics  of  the 
Koyna  region  have  not  been  extensively  studied  and, 
therefore,  are  only  partially  known.  The  Koyna  Reservoir  is 
located  a  few  kilometers  east  of  the  Western  Ghats,  a  1000m 
west-facing  escarpment  which  acts  as  the  continental  divide 
for  Peninsular  India  (Snow,  1974).  The  most  prominent 
^ Ibho logical  unit  in  the  area  is  the  Deccan  Trap  basalts 


Dube  obtained  his  crustal  model  from  the  body  wave  travel 
times  of  forty  Koyna  earthquakes  (mag  2-  4.0).  Although  the 
study  lacks  the  exactness  of  explosion  studies  like  those  of 
Kaila  et  al .  (1979)  and  Srivastava  et  al.  (1984),  it  is, 
however,  the  more  appropriate  model  for  this  study.  Since 
the  Dube  study  used  information  from  stations  distributed  in 
azimuth  in  the  Deccan  Trap  area  and  immediate  vicinity,  it 
gives  a  more  representative  indication  of  the  regional 
structure  than  the  others  studies  whose  information 
represented  only  a  few  localized  profiles  near  the  dam.  The 
Dube  model  is  also  supported  by  similar  results  from  the 
surface  wave  dispersion  study  of  Bhattacharya  (1981),  but 
Dube’s  model  was  selected  because  it  incorporated  the  work 
of  Tandon  (1973)  and  included  an  upper  layer  corresponding 
to  the  Deccan  Trap  basalts.  This  upper  layer  will  prove 
significant  in  affecting  the  short-period  waveforms  even 
though  it  is  very  thin.  It’s  significance  will  be  tested  by 
calculating  synthetics  for  a  three-layer  model,  model  B 
(Table  1),  without  the  upper  basaltic  trap  layer,  and 
comparing  them  to  the  synthetics  using  the  four-layer,  trap 
model . 

IV.  SYNTHETIC  SEISMOGRAMS 


THEORY 


To  investigate  the  effects  of  regional  structure  on  the 
propagation  of  seismic  waves  a  method  of  computing  synthetic 


seismograms  the  numbers  provide  a  convenient  way  of 
referring  to  specific  arrivals  whose  ray  paths  correspond  to 
those  in  figure  3. 

The  second  of  the  two  methods  to  be  used  is  a 
wavenumber  integration  technique.  This  approach  differs 
from  the  generalized  ray  theory  method  in  that  the 
calculations  are  done  in  the  frequency  domain  where  an 
integration  over  horizontal  wavenumber  or  slowness  is 


performed.  This  in  theory  allows  a  "complete”  solution  to 
be  calculated  resulting  in  a  synthetic  seismogram  with  all 
the  possible  ray  responses  such  as  multiple  reflections  and 
mode  conversions  included.  For  the  generalized  ray  method 
to  match  this  result,  an  infinite  number  of  rays  would  have 
to  be  specified.  The  theory  behind  the  wavenumber  technique 
is  described  in  Apsel  (1979)  and  Barker  (1984),  so  only  a 
brief  overview  of  the  method  will  be  presented  here. 

The  scalar  potential  wave  equation  in  cylindrical 
coordinates  is  solved  in  the  frequency  domain  using  the 
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where  ♦  ,  ♦  ,  and  X  ,  represent  the  P,  SV,  and  SH  potentials 
respectively.  A  represents  the  azimuthally  dependent  terms 


synthetics,  256  frequency  samples  were  computed  over  an 


interval  in  ray  parameter  from  .1  to  .4  sec/km.  These  were 


their  velocity  and  retrograde  motion.  These  waves  are 
designated  Rg  because  they  are  short-period  fundamental-mode 
Rayleigh  waves  which  travel  over  a  short  continental  path 
and  exhibit  normal  dispersion  (Bath, 1975).  On  the 
tangential  component  the  surface  waves  are  identified  as 
Love  waves  which,  in  this  case,  are  the  fundamental -mode 
transverse  waves  that  travel  in  the  Deccan  Trap  wave  guide 
at  velocities  approaching  the  S  velocity  (Bath, 1975). 

The  first  aspect  of  the  synthetic  plots  to  be  examined 
is  the  effect  of  varying  source  depth  on  the  generated 
waveforms  (Figures  4. 1-4. 5).  The  most  striking  observation 
is  the  rapid  decay  of  the  surface  wave  amplitude  as  the 
source  depth  increases.  At  a  source  depth  of  2  km  the 
surface  waves  are  the  largest  amplitude  waves  on  the  trace 
but  are  just  barely  identifiable  at  5  km  and  non-existent  at 
a  source  depth  of  10  km.  A  depth  effect  is  also  seen  in  the 
P  and  S  waves  where  the  duration  of  the  arrivals  get 
increasingly  longer  as  the  source  depth  increases.  As  the 
source  gets  deeper  the  time  differential  between  the  first 
and  last  major  arrival  gets  greater  causing  the  signal 
duration  to  lengthen.  The  reason  for  this  behavior  is  shown 
with  the  help  of  Figure  3.  The  last  major  arrival,  ray  #7, 
is  a  surface  reflection  and  a  source  depth  increase  would 
lengthen  its  travel  time  to  the  surface  and  therefore,  to 
the  receiver.  The  first  arrival,  on  the  other  hand,  changes 
from  the  direct  ray,  ray  #1,  to  a  higher  velocity  earlier 
head  wave  arrival  as  the  depth  increases. 


Another  thing  to  note  is  the  characteristic  ana 
distinctive  nature  of  the  major  S  wave  peaks  for  different 
source  depths.  These  major  S  wave  arrivals  are  identified 
by  numbers  in  Figures  4. 1-4. 5.  Specific  major  arrivals  can 
be  seen  moving  apart  as  the  depth  increases.  This  is 
particularly  apparent  in  rays  #4  and  #7  on  the  VSS  and  CLVD 
radial  traces  in  Figures  4.3  and  4.4  and  in  rays  #4,  #6,  and 
#7  on  the  VSS  and  VDS  tangential  traces  in  Figure  4.5.  The 
behavior  of  these  peaks,  that  are  easily  identifiable  with 
the  help  of  generalized  ray  theory  will  be  a  key  factor  in 
determining  the  depth  and  distance  of  the  real  events. 

One  last  thing  to  note  is  the  lack  of  any  systematic 
amplitude  change  in  the  body  waves  that  can  be  established 
as  a  function  of  depth.  This  is  explained  by  the  fact  that 
the  maximum  amplitude  depends  so  much  on  largest  S  wave  peak 
which  is  very  sensitive  to  angles  of  incidence  and 
constructive  and  destructive  interference  of  arrivals  caused 
by  structure  effects  at  this  distance. 

The  next  aspect  to  be  discussed  is  the  effect  of  event 
distance  on  the  synthetic  waveforms  which  is  most  useful  for 
determining  event  locations  (Figures  5. 1-5.5). 

Unfortunately,  only  estimates  on  event  distance  and  not 
location  can  be  made,  because  with  these  data  from  only  a 
single  station,  event  azimuth  cannot  be  accurately 
determined.  S-P  times,  the  time  difference  between  the 
initial  P  and  S  arrivals,  are  an  effective  way  of 
determining  event  distances.  The  S-P  times  as  seen  in 
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figures  5. 1-5. 5  range  from  14  seconds  at  110  km  to  18.5 
seconds  at  150  km.  The  initial  P  and  S  times  were  selected 
from  the  vertical  and  tangential  components  respectively. 

S-P  times  determined  in  a  similar  fashion  wJ 11  be  used  in 
estimating  the  distances  of  the  real  events. 

Another  significant  aspect  is  again,  the  variation  of 
the  S  wave  peak  amplitude.  One  cause  of  this  variation  is 
evident  in  the  S  wave  portion  of  the  vertical  and  tangential 
components  (Figure  6).  Here  the  interaction  between  ray  #6 
and  ray  #4  plays  an  important  role  in  shaping  the  S 
waveform.  As  the  distance  increases  the  order  of  the 
arrivals  changes  from  #6  arriving  before  #4,  at  110-120  km, 
to  simultaneous  arrivals,  at  130  km,  to  #4  arriving  before 
#6 ,  at  140-150  km.  This  wave  propagation  effect  of 
constructive  and  destructive  interference  significantly 
alters  the  relative  S  wave  amplitudes  and  character  of  the  S 
waves  on  these  components. 

A  different  type  of  distance  effect  can  also  be  seen  in 
the  S  wave  portion  of  the  seismograms.  The  timing  of 
post-critical  reflections  alters  the  nature  of  the  traces. 
Impressive  changes  in  arrival  #7,  the  surf ace-to-Moho 
reflection,  are  evident,  changing  from  a  small  pre-critical 
arrival  at  120  km  to  a  large  post-critical  reflection  at  a 
distance  of  130  km  (Figure  7). 

When  the  variation  in  source  mechanism  is  examined  it 
must  be  kept  in  mind  that  one  information  to  be  gained  is 
very  limited.  Due  to  the  fact  that  only  pure  VSS,  VDS, 


CLVD,  and  explosion  sources  are  calculated  at  a  single 
station  and  azimuth,  only  major  differences  will  be  of  any 
use  when  analyzing  the  real  events.  At  first  glance,  the 
only  characteristic  that  stands  out  is  the  fact  that  the 
amplitudes  of  the  radial  S  waves  are  smaller  than  the 
vertical  S  waves  for  VDS  mechanism  by  a  factor  of  thi^ee.  In 
comparison,  the  ratio  of  the  vertical  to  radial  S  waves  is 
about  one  for  the  VSS  and  CLVD  mechanisms.  The  variation  in 
source  radiation  pattern,  which  has  a  profound  affect  on  the 
seismic  waveforms,  causes  these  changes  in  amplitude. 

The  explosion  source  shown  in  Figures  4.2  and  4.4,  even 
though  not  directly  relevant  to  the  Koyna  earthquakes,  is 
still  useful  in  showing  wave  propagation  effects  for 
consider  iion  of  nuclear  discrimination  problems.  Since  the 
explosion  source  generates  no  S  wa/e,  it  is  also  useful  for 
examining  possible  P  to  S  conversions.  As  seen  in  Figures 
4.2  and  4.4,  P  to  S  conversions  (arriving  at  approximately 
the  same  time  as  the  S  arrivals  on  the  other  mechanisms)  do 
not  have  a  profound  effect  on  the  character  of  the 
seismogram.  Even  at  a  source  depth  of  2  km  the  major  P  to  S 
conversion  is  an  order  of  magnitude  smaller  than  the 
corresponding  S  waves  of  the  other  mechanisms.  The  same  is 
true  for  the  surface  waves. 

The  final  aspect  of  the  synthetic  seismograms  to  be 
examined  is  earth  structure  effects  which  are  model 
dependent.  To  do  this,  model  A  and  model  B  synthetics  are 
compared.  The  most  noticeable  difference  in  seismograms 
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calculated  from  model  B  (Table  1)  is  that  they  have  a  total 
lack  of  long-duration  surface  waves  even  at  a  source  depth 
of  2  km  (Figures  8. 1-8. 3).  This  effect  is  caused  by  the 
change  in  arrival  time  of  the  surface  waves  due  to  the 
existence  of  the  Trap  layer  in  model  A.  The  Rg  wave 
(Figures  8. 1-8.2)  travels  at  «•  3.0  km/sec  in  model  B  but  at 
only  *•  2.4  km/sec  in  model  A.  At  a  distance  of  120  km  this 
causes  the  model  B  Rg  wave  to  arrive  only  3  sec  after  the 
direct  S  while  the  model  A  Rg  wave  arrives  12  sec  after  the 
initial  S.  Thus,  model  B  surface  waves  interfere  with  S 
reflections  arriving  shortly  after  the  direct  S  arrival, 
while  model  A  surface  waves  are  clearly  visible  on  the 
trace,  arriving  long  after  the  major  S  reflections.  The 
same  is  true  of  the  Love  waves  in  Figure  8.3.  This  upper 
layer  also  has  the  effect  of  dispersing  the  surface  waves 
due  to  the  wave  guide  it  provides. 

Model  B  seismograms  are  less  complicated  in  terms  of 
general  appearance  than  those  of  model  A  because  there  is 
less  reverberation  in  the  structure  with  fewer  layers.  Also 
in  the  three-layer  model,  the  radial  P  waves  and  vertical  SV 
waves  are  large  while  the  opposite  is  true  for  the 
four-layer  model.  This  occurs  because  the  upper  trap  layer 
on  the  four-layer  model  acts  to  decrease  the  angle  of 
incidence  of  the  incoming  waves.  Thus,  the  P  waves  are 
larger  in  the  vertical  direction  and  the  SV  waves  are  larger 
in  the  radial  direction,  since  SV  motion  is  orthogonal  to 
the  P  wave  motion.  Another  difference  is  the  distance  where 
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arrivals  can  be  seen  moving  apart  as  the  depth  increases. 
This  effect  is  illustrated  in  Figure  9  where  the  real  events 
are  put  in  order  of  increasing  depth  based  on  the  synthetics, 
of  2,  5,  10,  and  15  km  source  depths,  assuming  these  phase 
identifications  are  correct.  The  arrival  times  of  both  rays 
#4  and  #7  were  corroborated  to  some  degree  by  the  vertical 
and  tangential  components  as  well. 

The  VSS  synthetics  are  shown  because  the  real  events 
seem  to  resemble  that  type  of  mechanism  more  than  any  of  the 
others  which  were  synthesized.  This  is  consistent  with  many 
of  the  mechanism  determinations  previously  reported  (Guha  et 
al.,1970;  Guha  et  al.,1974;  Langston, 1976 j  Rastogi  and 
Talwani , 1 980 ;  Langston, 1981 )  which  show  that  strike-slip 
faulting  is  common  in  the  area.  The  strike-slip  mechanism 
is  also  supported  by  the  amplitude  of  the  S  waves  seen  in 
the  real  events.  The  S  wave  amplitude  is  much  greater  than 
the  P  wave  amplitude  and  dominates  the  character  of  the 
seismogram.  This  was  seen  in  the  VSS,  but  not  in  the  VDS 
synthetics  (Figures  4. 1-4. 5). 

The  events  shown  in  Figure  9  seem  to  have  occurred  at 
distances  of  less  than  130  km  because  of  the  small  amplitude 
of  ray  #7  compared  to  the  large  amplitude  of  ray  #4 .  This 
amplitude  variation  was  discussed  in  the  synthetic  results 
section  and  is  shown  in  Figure  7  where  the  surface  to  Moho 
reflection  changes  from  of  small  amplitude  arrival  at  120  km 
to  a  large  amplitude  post-critical  reflection  at  130  km. 

The  S-P  times  of  the  real  events  also  indicate  distances  of 


less  than  130  km.  The  S-P  times  for  the  real  events  range 
from  14  to  17  seconds.  When  compared  to  the  synthetics  this 
indicates  event  distances  between  110  and  130  km. 

In  comparing  the  results  from  the  two  models  (Figures 
8. 1-8. 3),  it  is  determined  that  model  A  is  the  superior 
model,  thus  making  the  trap  layer  an  important  ingredient  in 
characterizing  regional  wave  propagation  of  the  Koyna  area. 
Many  differences  in  the  synthetics  generated  for  the  two 
models  make  it  clear  that  model  A  is  the  better  model.  The 
first  of  which  is  that  there  are  surface  waves  present  in 
some  of  the  real  events.  Their  long  duration  and  timing 
cannot  be  duplicated  in  the  three-layer  model,  even  at  a 
source  depth  of  2  km,  but  is  duplicated  by  the  four-layer 
model  with  the  addition  of  the  trap  layer.  The  longer 
duration  of  the  surface  waves  in  the  real  events  could  be  a 
result  of  dispersion  or  reflecting  and  scattering  of  waves 
by  lateral  heterogeneity  in  the  trap  layer. 

There  are  also  more  detailed  features  in  the  observed 
data  which  are  duplicated  more  effectively  with  model  A  than 
with  model  B.  The  surf ace-to-Moho  reflection  is 
post-critical  at  120  km  in  model  B  (Figure  8.1)  but  not  in 
model  A.  This  tends  to  support  model  A  because  the  real 
events  do  not  show  a  large  arrival  at  that  time  and  the 
events  are  in  that  distance  range.  The  amplitude  ratios  of 
P  to  SV  waves  on  the  vertical  and  radial  components  also 
agree  more  with  model  A.  The  P  wave  radial-to-vertical 
amplitude  ratio  is  much  greater  for  the  three-layer  model 


information  from  this  study  are  combined  with  information 
from  other  stations  in  the  area,  much  more  could  be  learned 
about  the  nature  of  the  event  source  parameters .  This  would 
help  increase  the  knowledge  of  the  geologic  and  tectonic 
relationships  of  these  reservoir  induced  events  from  which 
some  of  the  causes  could  eventually  be  determined. 
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TABLE  1 


MODEL  A 

Four-Layer  Model  Structure  Parameters 


*  *  *  *  *  **  *  *  *  *  *  *  5|c  *  *  *  He  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  s*c  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  * 


Vp 

(km/sec) 

Vs 

(km/sec) 

Density 
( S/cc ) 

Thickness 

(km) 

QPi 

QSi 

4.60 

2.65 

2.60 

1.2 

1000 

500 

5.78 

3.34 

2.70 

20.0 

1000 

500 

6 . 58 

3.80 

2.90 

18.8 

1000 

500 

8.19 

4.73 

3.20 

H.S. 

1000 

500 

******  *  S(C  *****  *  *******  *  *  *  *  *  **  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  * 


MODEL  B 

Three-Layer  Model  Structure  Parameters 


********* * * * ****** * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 


vP 

(km/sec ) 

Vs 

(km/sec) 

Density 

(g/cc) 

Thickness 

(km) 

QPi 

QSi 

5.78 

3.34 

2.70 

20.0 

1000 

500 

6 . 58 

3.80 

2.90 

18.8 

1000 

500 

8.19 

4.73 

3.20 

H.S. 

1000 

500 
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1  Attenuation  not  used  in  generalized  ray  calculations 
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Figure  2.3  Tangential  component  of  the  thirteen  digitized  Koyna 

events  of  December  1967  recorded  at  POO.  Major  phases 
are  denoted  by  S  and  L  (Love) .  Amplitudes  are  normalized 
and  given  in  microns. 


MODEL  A 


VERTICAL  (A  =  1  30KM) 


VERTICAL 


Synthetics  created  using  Model  A  showing  waveform  variatioi 
depth  of  radial  component  for  compensated  linear  vector  di 
(CLVD)  and  explosion  (EXPL)  source  mechanisms.  Note  the  S 
character  and  the  presence  of  surface  waves  at  2  km  depth. 
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lynthetics  created  using  Model  A  showing  waveform  variation  with 
listance  of  tangential  component  for  vertical  strike-slip  (VSS)  and 
ertical  dip-slip  (VDS)  source  mechanisms. 


the  wavenumber  synthetics  and  the  real  events  even 
though  only  seven  ray  paths,  as  shown  in  Figure  3, 
are  specified. 
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VERTICAL  RADIAL  TANGENTIAL 


RADIAL  (A  =  1  20KM 

li  2  68X10“* 
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MODEL  A  *  i2.5  sec  I  MODEL 


MODEL  A _ I  ,2  5  sec  I _ MODEL  B _ 

— Comparison  of  Modal  A  and  Model  B  tangential-component  synthetics 
(VSS  mechanism). 


Radial  component,  of  both  synthetic  and  real  events  as  a  function  of 
increasing  source  depth.  Depth  is  determined  by  the  separation  of 
S  and  sS  arrivals  (rays  *4  and  #7). 
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RADIATION  CHARACTERISTICS  OF  ELASTODYNAMIC  LINE  SOURCES  BURIED  IN  LAYERED 
MEDIA  WITH  PERIODIC  INTERFACES.  I  -  SH  WAVE  ANALYSIS 


Vijay  K.  Varadan+,  Akhlesh  Lakhtakia+,  Vasundara  V.  Varadan+  &  Charles  A.  Langston* 

^Laboratory  for  Electromagnetic  and  Acoustic  Research 
Department  of  Engineering  Science  and  Mechanics 
Department  of  Geosciences 
The  Pennsylvania  State  University 
UNIVERSITY  PARK,  PA  16802 


ABSTRACT 

Using  the  T-matrix,  or  the  extended  boundary  condition,  method  the  solution  for  a  class  of  problems 
involving  an  SH  line  source  in  an  elastic  waveguide  is  determined.  The  boundaries  of  the  waveguide  may  be 
periodically  corrugated  and  the  waveguide  may  be  embedded  between  elastic  media.  Numerical  results  ;xe  given 
for  a  seismically  interesting  case  of  wave  propagation  in  a  one-layer  crustal  model  over  a  mantle  hair-space  with 
a  corrugated  free  surface  representing  the  Basin  and  Range  topography  in  the  Western  United  States.  Analysis  of 
the  scattered  fields  at  the  surfxe,  and  of  the  fields  radiated  into  the  halfspace,  shows  complicated  field  behavior, 
even  with  sinusoidal  free  surface  corrugation.  These  results  are  directly  applicable  to  regional  wave  propagation 
and  scattering. 


118 


INTRODUCTION 

Scattering  is  known  to  be  a  general  characteristic  of  seismic  wave  propagation  and  undoubtedly,  it 
plays  a  iarge  part  in  the  observational  interpretations  of  Q  and  the  mechanism  of  anelastic  attenuation. 
Scattering  in  short-period  regional  phases,  such  as  Lg  and  Pg,  is  very  important  in  interpretations  of  the  seismic 
source  spectrum  and  radiation  pattern  discrimination  between  small  explosions  and  earthquakes.  For  example, 
transverse  L  motion  is  nearly  always  observed  from  explosion  sources,  yet  a  purely  isotropic  source  should 

o 

excite  far-field  motion  only  in  the  vertical  plane.  Either  the  source  is  not  isotropic,  which  would  imply  some 
common  form  of  tectonic  release,  or,  more  likely,  transverse  motion  is  due  to  the  scattering  by  the  geometric 
irregularities  of  the  crustal  waveguide.  Source  spectrum  estimates  from  these  phases  rely  on  the  assumption  of  a 
white  spectrum  for  propagation  effects  to  yield  interpretations  of  seismic  moment  and  comer  frequency,  but  as 
yet,  this  assumption  has  not  been  tested  for  validity;  scattering  must  play  a  role  in  shaping  the  seismic 
spectrum  of  these  phases. 

At  present  there  are  no  practical  numerical  techniques  to  compute  the  far-field  displacement  in  even 
simple  models  of  structures  containing  irregular  boundaries.  Although  techniques  which  incorporate 
plane-layered  elastic  structure  models  have  been  very  useful  in  routine  computations  and  data  interpretations, 
these  methods  are  inacapable  of  rigorously  addressing  even  the  simplest  of  laterally  heterogeneous  structures. 
Methods  are  needed  to  extend  these  computations  of  plane  layered  structure  to  include,  as  a  first  perturbation,  the 
effect  of  changing  boundary  geometry.  The  method  presented  in  this  paper  is  developed  for  geometries  of  two 
dimensions  and  SH  line  sources.  It  represents  a  necessary  first  step  in  dealing  with  waveguide  scattering  due  to 
heterogeneous  boundaries  and  will  be  extended  to  geophysically  interesting  point  source  problems  in  the  near 
future.  Nevertheless,  even  this  relatively  simple  problem  yields  insight  into  scattering  mechanisms  present  in 
the  crustal  waveguide. 

We  will  utilize  the  extended  boundary  condition  method,  also  known  as  the  T-matrix  method 
[Waterman,  1975],  to  investigate  the  radiation  characteristics  of  a  source  embedded  in  a  material  plate  bounded 
on  the  top  by  a  traction-free  surface,  and  on  the  bottom  by  a  homogeneous  half-space.  Both  surfaces  can  be 
taken  to  be  periodically  rough  in  the  theory  to  be  developed  here,  but  for  the  sake  of  simplicity,  we  shall  take 
the  lower  one  to  be  planar.  Discrete  mode  representation  of  the  fields  involved  will  be  used  to  solve  the  radiation 


problem.  Because  of  the  relative  simplicity  of  SH  wave  propagation,  the  features  of  this  technique  are  more 
clearly  displayed  without  the  complicating  effects  of  coupling  between  P  and  SV  wave  fields.  Hence,  in  this 
paper  we  are  interested  only  in  the  radiation  characteristics  of  a  SH-wave  line  source,  while  P-SV  line  sources 
are  treated  in  a  companion  paper  [Varadan  etal.,  1986]. 

Several  interesting  problems  can  immediately  be  examined  fiom  the  theoretical  treatment  presented. 
As  an  example,  both  interfaces  can  be  made  bimaterial,  thereby  allowing  one  to  consider  the  radiation 
characteristics  of  a  line  source  within  a  leaky  parallel -plate  waveguide.  Another  problem  of  interest  would  be  to 
have  both  interfaces  impermeable,  which  could  be  used  to  study  the  radiation  characteristics  of  a  line  source  in  a 
non-leaky  waveguide  of  varying  cross-section.  Yet  another  possibility  would  to  be  to  replace  the  line  source  by 
any  other  source,  provided  a  plane  wave  spectral  decomposition  of  the  field  radiated  by  it  is  possible  [Lakhtakia 
et  al.,  1985a],  Thus,  it  can  be  concluded  that  the  theory  presented  is  versatile  enough  to  be  of  use  in  solving 
problems  in  several  areas  of  research. 

THEORETICAL  DEVELOPMENT 

Consider  a  SH-wave  line  source  which  is  embedded  at  a  point  Tp  in  between  two  non-intersecting 
surfaces  Sa  and  Sb  whos*.  meui  planes  are  parallel  to  the  x  axis.  While  Sa  is  specified  by  the  function  y  =  h(x) 

>  yp  with  h(x)  possessing  a  finite  differentia]  coefficient  and  having  a  period  I  in  its  argument,  i.e.,  h(x±L)  - 
h(x);  for  the  sake  of  simplicity  we  consider  the  surface  Sb  to  be  the  plane  y  »  -  d  <  yp.  The  surface  Sa  is 
assumed  to  be  traction-free;  the  elastic  medium  in  between  the  surfaces  has  a  shear  modulus  of  |i  and  density  p, 
while  the  homogeneous  half-space  below  Sb  is  characterised  by  p.'  and  p'.  In  what  follows  the  ratio  p'/p  *  q,  k  = 
w[p/p] and  K  =  co[p'/p'l^;  an  exp[-iox]  harmonic  time  dependence  has  been  assumed;  i  and  j  represent  the  x- 
and  the  y -directed  unit  vectors,  respectively;  and  the  scalar  w  represents  z-directed  displacement  fields. 

On  applying  Huyghens'  principle  to  the  region  enclosed  between  Sa  and  Sb  in  the  manner  of 
Waterman  [1975]  and  Bostrom  [1980],  we  find  that  the  total  field 

w(r)  -  w‘(r)  +  (l/4i)  jSb  dsQ  nb(rQ)  •  [H^r  -  rj)  V()wb(r0)  -  w^  V^Odr  -  rj)] 

+  (L/4i)  JSads0  na(ro)  •  [wa(ro)  VoH0(k]r  -  rj)],  if  h(x)  >  y  >  -  d 
for  all  r  in  between  Sa  and  Sb,  while 


(la) 


with  Jn[*]  being  the  cylindrical  Bessel  function  of  order  n,  The  twin  sets  of  matrix  equations  (10a,b)  can  be 
solved  very  easily,  and  the  coefficients  A^a)  and  Bm(a)  computed  for  each  a. 

DISCUSSION  AND  NUMERICAL  RESULTS 

In  presenting  the  data  computed  using  the  present  method,  however,  we  shall  confine  ourselves  to  the 
surface  Sa  being  given  by  (12a)  and  Sb  being  the  plane  y  « -  d.  Throughout  this  section  now  we  shall  keep  a  =  d 
-  L,  p  -  2700  kg  m'3  and  \i  =  3.3075xl010  N  m'2;  the  line  source  is  placed  at  rp  =  (0,0).  Of  necessity,  the  sets 
of  equations  (10a,b)  have  to  be  terminated  at  some  |m|,  |n|  «  N,  it  being  made  sure  that  convergence  is  achieved 
at  some  value  of  N  for  each  a.  This  is  done  by  ensuring  that  the  unitarity  parameter  [Lakhtakia  et  at.,  1985b] 
for  each  a 

£W  <  N  I  Am(a)/4rt  Wa>|2  MYm(°)} 

converges  to  within  a  tolerance  limit  of  0.1%  as  N  increases.  Simultaneously,  it  is  also  ensured  that  Am(a)  and 
Bm(a)  converge  satisfactorily  as  well.  If,  in  addition,  a  sufficiently  large  number  (91,  here)  of  a 's  are  chosen  in 
the  range  |a|  5  kI L,  it  is  reasonable  to  expect  that  a  satisfactorily  converged  solution  for  the  total  problem  has 
been  obtained.  This  is  because  the  various  a 's  in  this  domain  do  not  interact  with  each  other  owing  to  the 
periodicity  of  the  problem.  It  may  be  mentioned  here  that  values  of  N  no  larger  than  13  turned  out  to  be 
adequate  for  the  calculations  presented  in  this  section. 

The  calculation  of  the  far-  zone  radiated  field  in  the  region  below  Sb  is  quite  straightforward.  If,  for  a 
given  m  and  a,  Ym(a)  be  real  positive,  then  at  an  angle  9  -  tan'^ct/Yj^fa)}  with  respect  to  the  y-axis,  the 
radiated  field 

w2(r,0)  {A^cOAlrriYjj/a)}  exrriKm“(a)  •  r]  as  Kr  ->  «*.  (14) 

Shown  in  Figs.  1  and  2  are  these  radiation  amplitudes,  j  Am(a)/4r,iYm(a)|  as  functions  of  the  angle  9  - 
sin'1(a+2rvm/L ),  |a+2~m/L  |  £  K,  |m|  <■»,  for  roughness  b/L  -  0.10  at  the  normalised  frequencies  kL  =  5.0 
and  10.0.  The  medium  below  Sb  has  a  density  p'  -  3200  kg  m'3  and  a  shear  modulus  ja'  -  6.48xl010  N  in*2. 
The  chosen  values  for  the  densities  and  the  shear  moduli  correspond  to  shear  wave  velocities  of  3500  and  4500 
m  s’1,  appropriate  for  a  one-layer  continental  crust  model.  Normalized  frequencies  used  correspond  to  shear 
waves  of  periods  10  and  5  s,  respectively.  The  corrugation  period,  L,  is  equal  to  the  mean  crustal  thickness,  and 


the  choices  for  the  corrugation  amplitude  are  based  on  the  average  wavelength  of  the  topography  encountered  in 
the  Basin  and  Range  province  of  the  western  United  States.  Thus,  b/L  =  0.1  corresponds  to  topography  with  3 
km  amplitude  and  30  km  wavelength.  The  examples  show  wave  propagation  for  waves  with  wavelengths  shorter 
than  the  topography  and  appropriate  for  intermediate  period  seismic  data. 

Two  points  need  to  be  noted  here.  Firstly,  the  use  of  exp[i(a  +  2rm/L)x]  in  (9)  instead  of  exp[ikn*(a) 
•  rj  has  allowed  us  to  compute  for  roughnesses  b/L  >  0.072,  the  classic  Rayleigh  limit  [Miliar,  1973;  Lakhtakia 
et  at..,  1985c]  for  surface  scattering  problems.  This  has  been  established  by  several  authors  [e.g.,  Masel  et  al., 
1975,  Goodman,  1977;  Cain  et  al.,  1986]  to  whom  the  interested  reader  is  referred.  Secondly,  the  computed 
amplitudes  are  somewhat  asymetricai  about  9  »  0°.  This  is  because  the  upper  surface  is  asymmetrical  about  the 
location  of  the  line  source,  due  to  the  former's  specification  as  y  «  a  *  fc  ji>.(2jix/L). 

An  important  characteristic  of  Figs.  1  and  2  is  the  presence  of  sharp  extrema  in  the  radiation  pattern 
plots.  Not  counting  the  traction-free  surface  Sa,  for  the  moment,  these  extrema  can  occur  for  at  least  two 
reasons.  Firstly,  the  presence  of  the  Brewster  anomaly  [Weeks,  1964],  whenever 

a  +  2rtn/L  -  K{[q2  -  (k/K)2]/[q2  -  1]}*1/2  ,  |n|  -  0,1, 2, 7  (15) 

for  a  given  a  and  n,  is  such  that  Rn(a)  is  real  positive,  maximizes  transmission  across  a  planar  bimaterial 
interface.  Also,  if  a  +  2run/L  -  ±K,  then  a  reflectionless  anomaly  [Cain  et  al.,  1986]  appears  for  a  planar 
bimaterial  interface. 

In  order  to  consider  yet  another  reason  for  the  extrema  in  these  figures,  let  us  consider,  first,  a  planar 
Sa.  When  the  .wo  waves  exp[ikn±(a)  •  (r  -  rp )]  launched  by  the  line  source  first  arrive  at  S^,  their  phases  differ 
by  2aBn(a).  It  2aBn(a)  is  equal  to  an  odd  multiple  of  ±Jt  they  interefere  destructively,  but  if  2aBn(a)  is  an  even 
multiple  of  irt  then  constructive  intereference  takes  place.  Consequently,  minima  and  maxima  may  occur.  For 
example,  when  b/L  -  0  in  Fig.  2,  a  minimum  is  recorded  in  the  radiation  pattern  when  0  -  25.5H,  a/K  * 
0.4297745  and  2aBQ(a)  -  3tt.  When  the  upper  surface  becomes  periodically  rough,  such  simple  arguments  do 
not  explain  all  the  possible  extrema  [Lakhtakia  et  al.,  1985a).  It  may  be  surmised,  however,  that  the  mutual 
intereference  of  the  n  *  0  modes  in  (I0a,b)  makes  the  field  in  the  region  between  Sa  and  very  complicated, 
giving  rise  to  an  extremely  complex  radiation  pattern.  Due  to  periodicity,  nevertheless,  a  very  clearcut  extremum 
occurs  whenever  a  +  2rm/L  equals  either  ±k  or  ±K.  Such  an  extremum  is  called  the  Rayleigh-Wood  anomaly 


and  has  been  discussed  elsewhere  [Waterman,  1975;  Goodman,  1977;  Lakhtakia  et  al.,  1985a]  in  greatei  deuir 

In  Fig.  3  is  shown  the  radiation  pattern  when  (A)  q  *  1.9591837  and  K  -  (7/9)k,  and  (B)  q  -  1  and  k  = 

K.  Both  plots  are  at  a  normalised  frequency  kL  -  15  when  b/L  -  0.15.  To  be  noted  is  the  significant  complexity 
of  the  plot  in  Fig.  3(A)  when  compared  with  that  in  Fig.  3(B).  In  the  latter  case  the  media  on  either  side  of  the 
surface  Sb  are  identical,  which  means  that  the  extrema  due  to  the  back-and-forth  bouncing  of  the  waves  in 
between  Sa  and  Sb  are  no  longer  there. 

An  important  quantity  measured  for  seismological  applications  is  the  field  wa  induced  on  Sa.  Shown 
in  Fig.  4  are  the  computed  values  of  the  magnitude  of  wa  at  locations  x/L  *  ±1,  ±10  and  ±50  as  the  normalized 
frequency  increases  from  kL  -  10  to  kL  -  30  in  steps  of  0.2;  the  other  parameters  are  b/L  -  0.10,  q  -  1.9591837 
and  K  -  (7/9)k.  Equation  (9)  was  used  to  calculate  wa  with  a  3 1  -  point  Gauss-Legendre  quadrature  implemented 
for  the  integral  over  a.  The  asymmetry  in  these  computations,  again,  is  due  to  the  fact  that  Sa  is  asymmetric 
with  respect  to  the  line  source.  Although  the  plots  are  very  complicated,  generally  speaking  it  appears  that  wa 
decreases  somewhat  as  kL  increases.  Furthermore,  there  is  some  indication  of  a  trend  for  wa  to  decrease  as  the 
location  |x|  on  Sa  is  increased. 

In  conclusion,  the  examples  shown  in  this  paper  have  been  oriented  towards  intermediate  and 
short-period  seismic  wave  propagation  in  contendnental  crustal  structures.  Application  to  seismic  observations 
and  calculation  of  synthetic  seismograms  is  deferred  for  later  work.  The  important  result  is  that  this  technique 
makes  such  studies  feasible.  The  method  is  most  efficient  for  long  wavelength  interactions,  but  is  not 
practically  restricted  until  very  high  normalized  frequencies  are  attained.  The  limitation  on  interface  corrugation 
is  topography,  and  not  the  Rayle'.gh  ansatz  of  representing  each  surface  field  by  Bloch-Fourier  expansions. 
Although  a  SH  line  source  has  been  considered,  the  formalism  is  amenable  to  P-SV  line  sources  [Varadan  et  al., 
1986]  and  point  sources  as  we’1 
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Fig.l  Normalized  field  amplitude  as  a  function  of  0,  the  angular  displacement  from  the  y-axis,  in  the  medium 
below  Sb  due  to  a  line  source  located  at  rp  -  (0,0).  The  surface  Sa  is  y  -  L{  1  +  (0.10)  sin(27tx/L)},  while  the 
I  surface  Sb  is  the  plane  y  =  -  L.  The  relevant  parameters  are  kL  -  5.0,  q  -  1.9591837,  and  K  -  (7/9)k.  For 

comparison,  the  dashed  line  shows  the  same  calculations  for  b/L  -  0.  To  get  the  actual  amplitude,  multiply  the 
values  read  from  the  graph  by  0.126033  for  b/L  =  0,  and  by  0.102563  for  b/L  -  0.10. 
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Fig.4  Magnitude  of  wa  computed  at  locations  x  on  Sa,  due  to  a  line  source  located  at  Tp  -  (0,' 
of  the  normalized  frequency  IcL.  The  surface  Sa  is  y  -  L{1  +  0.10  sin(2jrx/L)},  while  the  si 
plane  y  -  -  L.  The  other  relevant  parameters  are  q  -  1.9591837,  and  K  -  (7/9)k.  (A)  x  -  +  L 
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ABSTRACT 

A  method  for  determining  the  elastodynamic  (P  and  SV  waves)  radiation  characteristics  of  finite-size 
sources  buried  in  horizontally  layered  media,  having  periodically  corrugated  interfaces,  is  described.  In  particular, 
the  example  problem  chosen  to  illustrate  the  procedure  is  as  follows:  a  solid  plate  lies  on  top  of  an  solid 
half-space;  the  solid-solid  interface  has  been  taken  to  be  planar,  but  traction-free  conditions  prevail  on  the  other 
boundary  of  the  elastic  plate,  which  surface  is  also  periodically  corrugated;  and  the  source  has  been  taken  to  be 
an  isotropic,  P-wave  line  source  located  inside  the  elastic  plate.  The  technique  presented  utilizes  the  planewave 
spectral  decomposition  of  the  relevant  fields  within  the  framework  of  the  extended  boundary  condition,  or  the 
T-matrix,  method.  Since  the  T-matrix  method  is  a  matrix  approach,  it  is  very  attractive  computationally,  and  is 
certainly  more  tractable  than  a  method  based  on  the  direct  solution  of  the  integral  equations  involved  in  the 
scattering  problem.  Numerical  results  are  given  to  delineate  the  various  features  of  the  field  diffracted  into  the 
elastic  half-space,  as  well  as  those  of  the  displacement  field  induced  on  the  traction-free  boundary  of  the  elastic 
plate.  The  specific  example  examined  is  directly  related  to  regional  wave  propagation  in  a  continental  crustal 
waveguide. 


INTRODUCTION 


In  the  companion  paper  [Varadan  et  al.,  1986],  a  method  was  developed  to  treat  the  scattering  of  the 
radiated  fields  of  a  SH-wave  line  source  embedded  in  a  two-dimensional  elastic  waveguide  bounded  by  corrugated 
interfaces,  which  can  lie  either  between  elastic  media,  or  have  homogeneous  boundary  conditions  prescribed  on 
either  or  both  of  the  interfaces.  Problems  of  this  type  are  important  in  various  fields  such  as,  noise  control, 
meteorology,  as  well  as  seismology. 

In  this  paper,  we  utilize  the  extended  boundary  condition,  or  the  T-matrix,  method  [Waterman,  1975]  to 
treat  the  problem  of  P-  and  SV-  wave  propagation  in  a  waveguide  excited  by  a  line  source.  This  is  a  necessary 
extension  of  the  Technique  to  treat  problems  involving  seismic  wave  p  opagation  at  regional  and  teleseismic 
distances  where  the  effect  of  waveguide  imperfections  is  of  great  interest.  A  motivation  for  obtaining  this 
solution  is  to  understand  crustal  waveguide  scattering  of  regional  waves,  such  as  L  and  P  excited  by  nuclear 
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explosions. 

A  simple  scattering  arrangement  is  chosen:  a  solid  plate  lies  on  top  of  an  elastic  (solid)  half-space.  The 
solid-solid  interface  has  been  taken  to  be  planar,  but  traction-free  conditions  prevail  on  the  other  boundary  of  the 
elastic  plate,  which  surface  is  also  periodically  corrugated;  and  the  source  has  been  taken  to  be  an  isotropic,  P 
wave  line  source  located  inside  the  elastic  plate.  Tit?  parameters  and  dimensions  of  the  model  chosen  hf.se  are 
motivated  by  the  geophysically  interesting  situation  of  a  continental  crust  over  a  mantle  half-space.  The 
corrugations  of  the  free  surface  simulate  the  topography  of  the  Basin  and  Range  province  of  the  western  United 
States,  as  has  been  alluded  to  in  the  companion  paper  [Varadan  et  al.,  1986],  It  should  be  noted  that  such  a 
simple  scattering  geometry  merely  exemplifies  more  complicated  arrangements:  there  can  be  several  elastic 
(solid  or  fluid)  plates  atop  each  other;  some  or  all  of  the  various  interfaces  can  be  periodically  rough;  and  the 
choice  of  the  source  can  be  purely  arbitrary,  provide  the  field  radiated  by  it  can  be  expressed  [Clemmow,  1956] 
in  terms  of  a  planewave  spectrum,  discrete  or  continuous. 

The  method  of  solution  proposed  here  is  not  an  extension  of  the  one  used  elsewhere  by  us  [Lakhtakia  et  a!., 
1985a],  where  all  relevant  fields  were  expanded  in  terms  of  a  discrete  planewave  spectral  representation,  the 
T-matrix  of  each  interface  was  separately  computed,  and  the  inter-surface  distance  was  taken  into  account  using 
co-ordinate  transfer  [Cain  et  al.,  1986]  or  planewave  propagation  matrices.  Instead,  here  the  extinction  theorem 
of  the  extended  boundary  condition  method  [Waterman,  1975;  Barber  and  Yeh,  1975]  is  applied  to  elastic  plate 
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and  the  relevant  fields  are  decomposed  into  a  continuous  planewave  spectrum.  The  resulting  theoretical  treatment 
is  then  converted  into  a  set  of  matix  equations,  which  can  then  be  solved  on  a  digital  computer.  This  process 
gives  rise  to  a  more  compact  formulation  as  compared  to  that  which  could  be  obtained  using  the  method  of 
Lalchtakia  et  al.  [1985a].  Numerical  results  are  given  to  delineate  the  various  features  of  the  field  diffracted  into 
the  elastic  half-space,  as  well  as  those  of  the  displacement  field  induced  on  the  traction-free  boundary  of  the 
elastic  plate. 

PRELIMINARIES  AND  NOTATION 

Shown  in  Fig.  1  is  a  schematic  of  the  problem  to  be  solved  here.  The  stress-free  boundary  condition 
prevails  on  Ta,  which  surface  is  determined  by  the  relationship  y  -  h(x)  >0,  <  x  <  where  h  is  single 

valued  and  has  a  period  L  m  its  argument,  i.e.,  h(x)  -  h(x  +  L).  On  the  other  hand,  the  surface  Tb  is  taken  to  be 
the  planar  interface  y  -  -  d  <  0,  -«><*<<»,  and  the  two  surfaces,  tfcv .  may  never  meet.  Interposed  in  between 
the  two  surfaces  is  a  homogeneous,  isotropic,  perfectly  elastic  solid  of  density  p  and  Lame  constants  X  and  p, 
while  the  semi-infinite  region  below  Tb  is  occupied  by  a  similar  elastic  solid  whose  constants  are  p',  X  and  p . 
Furthermore,  a  source  is  embedded  at  the  location  r^,  chosen  such  that  h(x)  >  yp  >  -  d,  Vx.  It  is  the  aim  of  this 
paper  to  obtain  the  field  present  in  the  primed  medium  below  Tb  due  to  this  source. 

Two-dimensional  wave  propagation  in  elastic  media  is  characterized  by  two  types  of  wave  motion:  the 
longitudinal  or  P  wave  and  the  vertically-polarized  shear  or  SV  wave.  In  general  mode  conversion  between  P 
and  SV  waves  takes  place  at  boundaries;  hence,  it  is  necessary  to  consider  the  two  wave  types  simultaneously. 
In  what  follows  a  combination  of  vector  and  tensor  notation  will  be  used,  and  which  will  be  developed  along 
with  the  theory.  The  position  vector  r  can  be  interchangeably  represented  as  either  (x,  y)  or  (rj,  rj).  Greek 
subscripts  (except  k)  are  to  be  only  assigned  values  of  1  and  2,  while  the  summation  convention  holds  for 
repeated  Greek  subscripts.  An  exponential  time  dependence  exp[-icotj  is  implicit  in  this  treatment 

Particularly  with  reftience  to  the  unprimed  medium,  the  displacement  and  the  stress  Green's  tensors, 
denoted  by  G(r  -  an  dl(r-  r^),  respectively,  together  sarisfy  the  differential  relations  [Achenbach,  1973): 

^"^~oJ3y  +  ^cc y  “ '  *  ro^’  ^ 

*^"ol{3y  +  P^(J^aY  +  ^a^f3Y-^  ^ 

with  both  G  and  I  also  satisfying  suitable  radiation  conditions  [Tan,  1975]  as  well  as  the  Floquet  conditions 


(2a) 


[Lakhtakia  et  al .t  1984b]  due  to  the  periodicity  of  Ta.  These  Green’s  tensors  are  represented  as: 

GaY<r-ro>“*LI^dK  DKGor/K-r-ro)' 

(r  ‘  ro>  “  ;VL^dK  Dk  T  r-  ro>‘ 

where 

1/Dk  -  1  +  5(k  -  jc/L)  +  8(k  +  tc/L),  p>r ) 

and  the  basic  variable  k  varies  from  -jc/L  to  nJ L  along  the  real  line,  and  8(k  -  k.')  being  the  Dirac  det'a.  Fv  :he 
solid  medium  with  density  p  and  Lame  constants  \  and  (i,  pertinent  to  the  present  problem,  the  foilowrr? 


quantities  are  now  defined: 

kxi~k  +  2jtn/L,  n  €  [0,  ±1,  ±2,  ±3,  ...],  (3a) 

*n  ■  +^2  *  Si2!:  CTn  "  +^2  -  (3b) 

Pn±  "  XKn  ±  y*n  ~  ^ln*  +  >p2n *  Sn*  "  XKn  ±  > °n  *  xS  In*  +  >  S2n*  (3c) 

kp  -  cVCp  ;  cp  «  [(X  +  2|i)/p] 1/2  (3d 

ks*:a/cs:cs“^yPl1/2-  (30 


Associated  with  each  k  are  also  the  various  free-space  Green’s  functions,  G^k,  r,  and  T^tc,  r„  yQ). 
The  free-space  displacement  Green's  function  can  be  most  easily  set  down  as  [Lakhtakia  et  al.,  i  84b]- 
Gcr/K* r-  ro)  "  Zn  (&s2/2pco2anL)  expfiS,*  •  (r  -  rQ)]  + 

+  Zn  (ikp2/2PO)2rtnL)  Uan1^  Uyn1*  “Pt^n*  *  <r  *  ro>^  <4' 

where 

^an*5  “  U^p)  Pan~*  ^anS_  ”  eaP  ^pn~’  (5a,b) 

Ea{3  *s  311  urdt  Asymmetric  tensor,  and  the  superscript  V  (resp.  '-')  is  to  be  used  if  y  >  yQ  (resp.  y  <  yQ).. 
Corresponding  to  this  G^x,  r,  is  the  stress  Green's  function  given  by  [Lakhtakia  et  al. ,  1984b] 

Totf*  r-  ro>  -  Zn  (^W^U  %S±  exp[iSn±  •  (r  -  r^]  + 

4  Sn  (ikp2/2pco27tnL)  U^P±  T^P*  exp[iPn±  •  (r  -  r^],  (6) 

where 

TaPaP±  -  ^  (wa3  +  ^  ucmP±  L'»nP±l-  (7a) 

W*  -  V  l  ^  U|)as±  ♦  Sp*  llma ).  (7b) 

It  should  be  noted  that  k  is  implicit  in  the  definitions  (3),  (5)  and  (87;  moreover,  by  putting  primes  on  the 
various  quantities  defined  above  (except  k),  those  quantities  can  be  made  to  refer  to  the  primed  medium  below 


APPLICATION  OF  THE  EXTINCTION  THEOREM 
The  application  of  Huyghens’  principle  to  the  field 
leads  to  [Waterman,  1975;  Varadan  et  al.,  1986]: 

[T •  ro>  uba(ro) 
+  /ra  * o  [Tofl^r-ro)uaa(ro> 


uaa(r)“-7C/Jn/LdK  uaa(K’r)’  (1  lb) 

etc.  Finally,  the  limiting  principle  [Waterman,  1971;  Waterman,  1975;  Barber  and  Yeh,  1975]  of  the  T-matrix 
method  has  to  be  employed.  Expansions  (4)  and  (6)  are  substituted  in  (1  la)  and  inner  products  of  (11a)  are  taken, 
with  exp[-iPn~(K)  •  r]  and  expf-uj^ic)  •  r] ,  over  the  domain  |x|  £  L/2  on  both  the  surfaces  Ta+  and  (see 
Fig.  1).  Of  the  four  inner  products  thus  taken,  those  over  are  straightforward  and  need  no  explanation.  But 
the  reason  for  taking  the  inner  products  over  ra+,  rather  than  on  Ta,  can  best  be  explained  by  the  concept  of 
analytic  continuation  put  forth  by  Waterman  [1971;  1975]  and  to  which  the  interested  reader  is  referred. 

As  a  result  of  the  foregoing  manipulations,  the  twin  relations  of  (11a)  simplify  to  the  four  equations, 
exp[-lPrf(K)  *  rp)  “  Cp  /ra  ^o  TopnP±  “aa^o*  nap(ro>  expI-i^K)  *  rj 

'  Cp  /rb  ^o  frctpn^  uba(K-ro)  +  Uctn1*  ^ap^o*]  "bp^  expt-iP^fK)  •  r^; 

Vk  ;  |k|  ^  ti/L,  Vn  :  |nj  ^  oo,  (12a,b) 

and 

0  =  Cs  /ra  ^o  TocPnS±  “aa^o*  nap(ro>  Mpf-iS^*)  *  ^ 

'  f  r~b  ^o  t^aPn  ~  “ba^^  +  ^ccn5*  TbaP^lc,rcP-l  nbP^ro^  axPt‘^n~(K^  *  ro^: 

Vk  ;  |k|  ^  rt/L,  Vn  :  |n|  £  °°,  (12c,d) 

where 

Cp =  i71  kp2/2pco2L,  C$  =  itt  ks2/2pco2L.  (13a,b) 

These  four  equaccas  suffice  to  solve  for  the  fields  induced  on  the  two  boundaries,  and  are  a  direct  consequence  of 
the  application  of  the  extinction  theorem.  Their  solution  can  be  attempted  using  the  integral  equation  methods 
[sg-»  Fokkema,  1979]  but,  instead,  here  the  four  equations  grouped  together  as  (12)  will  be  converted  into  a 
system  of  matrix  equations  after  using  the  boundary  conditions  on 

THE  MATRIX  EQUATIONS 
In  (12)  the  boundary  conditions 

Tbap(K,r)  nbpW  “  x'cxp(K’r>  nbp<r)  >  uba(K>r>  "  “'a**’1*  r  6  rb>  d4) 

appropriate  for  the  solid-solid  interface  should  be  substituted,  with  u'  and  x'  being  the  displacement  vector  and 

the  stress  tensor,  respectively,  of  the  field  radiated  into  the  region  below  The  plane  wave  spectral 
decomposition  of  these  two  quantities  is  given  by 


Further  simplification  by  p 
Vic  :  |k|  £  tt/L,  \ 
expt-iPjfOO  •  rp]  -CpZ 


(19atb) 


and 

0  “  Cs2m  {BmP  TaPnS±  UamP+  “pnm^n-  V> 

+  {Bm  ^apn1-  ^amS+  ^Pnm'^n’  cm^ 

*  LCs  \P  [Ta2nS±  u’o m1’'  +  T oSn^  Uom^W'^n  +  *'n)dl 
-  LCs  V  tTa2nS±  u'anS'  +  T a2nS'  Uan^^P^n  +  aY,)dl-  09c, d) 

where  the  integrals, 

N  lnm±(cl'  Q)  "  (2jt/L)  Kn-m)/(±q  -  Q))  X  2nm±(q.  Q).  (20a) 

and 

N2nm±(c!'  Q)  *  -Ul^  dt  exp[-i2(n-m)jnt/L]  exp[-i(±q  -  Q)h(x)].  (20b) 

If  the  traction-free  surface  Ta  be  sinusoidal,  i.e., 

h(x)  -  a  -  b  cos(2jtx/L),  (2ia) 

then  the  integral  (20b)  can  be  evaluated  in  terms  of  the  cylindrical  Bessel  functions  Jn(*)  as 

^nnfta-Q)  *  *,D  m|  L  expt-i(±q  -  Q)aj  J^fr+q  -  Q]b),  (21b) 

w‘th  N  lnm±^'  0)  still  defined  as  in  (20a);  the  only  special  case  arises  for  q  -  Q,  when 

N2nm+(<1-  <l)  *L8nm:  N  lnm+^  <D  »  ’i7lb  $n,m+l  *  5n.m-l]>  (21c,d) 

5am  being  the  Kronecker  delta. 

DISCUSSION  AND  NUMERICAL  RESULTS 

Of  necessity,  the  system  of  equations  in  (19)  has  to  be  truncated  at  some  |m|,  |n|  *  N,  for  each  k,  care  being 
taken  however  to  ensure  that  convergence  is  achieved  at  some  value  of  N.  This  is  done  by  ensuring  that  the 
unitarity  parameter  [Lakhtakia  et  al.,  1985b]  for  each  k, 

Zm  (*■' +  2H')  |ATnP(»c)I2  ^e{re’m(K)}+  |A'|Ams(iO|2 Re{a’m(K)}  (22) 

converges  to  within  a  preset  tolerance  of  ±0.1%  as  N  is  increased.  Simultaneously,  the  various  field  coefficients 
AjjjP’5  and  BmP's  should  also  converge  within  similar  error  bounds. 

It  should  be  noted  that  if  rcm(tc)  for  a  given  m  and  k  be  real  positive  then  at  an  angle  8  -  tan'1{K/nm(K)} 
with  respect  to  the  y-axis,  the  P  wave  component  of  the  radiated  field  in  the  medium  below  iy 
u'a(r)  AmP(K)U’amP'(,c)  exp[-iFm'(<)  •  r]  as  k'y  -» 


(23a) 


JKj? 


for  several  test  cases  were  carried  out.  In  Fig.  4,  the  radiation  amplitudes  at  kpL  -  7.5  are  plotted  for  (a+d)/L  - 
2.0,  b/L  -  0.05,  c'p  -  4500  ms‘',c's  -  8000  m  s"1  and  p'  -  3200  kg  m‘^;  the  separation  between  the  line 
source  and  assumes  two  values:  (a)  d/L  ■  1.0,  and  (b)  d/L  ■  1.5.  In  Fig.  5  the  calculations  of  Fig.  4  are 
repeated,  but  for  c'p  ■  Cp,  c’s  -  cs  and  p'  *  p.  Finally,  the  computations  of  Figs.  4  and  5  were  re-done  for  b/L  « 
0.0,  and  the  resulting  radiation  patterns  are  shown  in  Figs.  6  and  7,  respectively.  Thus,  all  relevant 
combinations  of  the  material  properties  and  the  geometry  of  r*a  have  been  used  in  Figs.  4  -  7  to  enable  one  to 
explain  the  features  mentioned  above. 

With  reference  to  Fig.  7,  First,  the  simplest  of  these  four  sets  of  calculations  (rtn  -  rc'n,  etc.),  consider  a  P 
wave,  with  a  wave  vector  Pn+(K),  launched  by  the  source  and  which  hits  the  plane  surface  Ta.  As  a  result,  a  SV 
wave  with  the  wave  vector  8^  Sm*(ic)  is  excited.  Because  of  Fresnel's  laws  [Lakhtakia  et  d.,  1984b;  Lakhtakia 
et  al,  1985a],  this  SV  wave  propagates  downward  at  an  angle  0,  where 

|9|  -  sin*1  {cs  |Pln+(K)|}  (24) 

with  respect  to  the  -y  axis.  If,  however,  |Pjn+(K)|  £  kp,  then  that  P  wave  carries  no  energy  in  the  +y  direction: 
consequently,  the  SV  waves  propagating  downwards  at  angles  |0|  >  sin-1  {cs/cp}  -  35.7  deg  may  not  transport 
any  energy  in  the  -y  direction.  It  is  conjectured,  therefore,  that  is  the  reason  why  the  SV  wave  amplitudes  in 
Fig.  7  quickly  become  negligible  around  this  limiting  angle.  ( In  Figs.  3  and  6,  the  corresponding  limit  will  be 
around  |0|  »  sin  *  {c's/Cp}  «  48.6  deg.)  When,  however,  Ta  is  periodically  corrugated,  as  in  Figs.  2  and  5,  the 
launched  P  wave  with  the  wave  vector  Pn+(K)  excites  all  SV  waves  with  wave  vectors  Sm*(K),  and  S'm*(K), 
thereby  extending  the  angular  spread  of  the  SV  wave  radiation  pattern. 

Considering  Fig.  7,  again,  it  is  known  [Achenbach,  1973]  that  a  Rayleigh  surface  wave  can  propagate  on 
the  traction-free  Fa  with  a  x-directed  wavenumber  ±kR  given  by 

kR  -  kj-ll+vMO.862  +  1.14V]*1;  0  5  v£  0.5,  (25a) 

where  v,  the  Poisson's  ratio,  is  defined  as 

2v  -  [(Cp/cs)2  -  2M(Cp/cs)2  - 1]*1.  (25b) 

Again,  if  attention  is  fixed  upon  any  k  :  |k]  <,  rc/L,  for  which  a  certain  k^k)  equals  ±kR,  then  the  excitation  of 
the  surface  wave  on  Ta  will  exhibit  itself  in  the  radiation  patterns  as  an  anomaly.  This  is  what  happens  for  ±0  - 
37.069  deg,  when  K+jfkj  sin©)  equals  ±kR,  and  the  spikes  in  the  SV  radiation  patterns  in  Fig.  7  exhibit  this 
phenomenon.  Using  the  methods  outlined  by  Lakhtakia  et  d.  [1984b;  1985a],  several  of  the  anomalies  in  Fgis. 


2  •  6  may  also  be  shown  to  be  attributed  to  the  excitation  of  the  Rayleigh  surface  waves  on  Ta. 

Finally,  and  again  with  reference  to  Fig.  7,  it  is  seen  that  the  P-wave  radiation  pattern  records  several 
broadened  extrema.  The  reason  for  these  is  as  follows:  The  P  wave  with  spatial  dependence  exp[iPn+(K)*(r  •  fp)] 
hits  Ta  and  is  reflected  back,  partly  ar,  a  P  wave,  with  a  change  of  phase,  say 
{^n(K)*exp[iPn'(K)*r]*exp[-iPn+(K)Tp]*exp[i2a7tn(K)]};  it  then  interferes  with  the  P  wave  having  the  spatial 
dependence  exp[iPn"(K)*(r  -  Tp)]  which  also  has  been  launched  by  the  line  source.  The  extrema  of  the  function 
{^(k)  exp[-iPn+(K)Tp]  expfLTa^fK)]  +  exp[-iPn'(tc)Tp]},  therefore,  cause  the  broadened  extrema  of  the  P-wave 
radiation  patterns  in  Fig.  7.  When  Ta  is  periodically  corrugated,  as  in  Fig.  5,  these  extrema  become  more 
numerous  and  sharper,  since  each  P  wave  incident  on  Ta  is  reflected  back,  partly,  as  an  ensemble  of  P  waves 
with  spatial  dependences  {^(k)  expfiPm*(K)*r]  exp[-iPn+(K)Tp]  exp[iann(K)]  exp[iartm(K)].  A  further  increase 
in  the  numbers  of  such  extrema  comes  in  if  is  a  bimaterial  interface,  as  in  Figs.  2  -  4,  and  6,  when  similar 
arguments  hold  for  the  multiply  reflected  P  and  SV  waves  inside  the  unprimed  medium.  In  addition,  because  of 
the  involvement  of  factors  like  exp[iartn(K)]  and  expOdji^K)],  etc.,  the  dependence  of  these  extrema  is  sensitive 
to  the  location  yp  of  the  line  source,  as  may  be  observed  by  comparing  'he  two  P-wave  radiation  patterns  in 
Fig.  7. 

Comparison  of  Figs.  5  and  7  reveals  yet  another  class  of  anomalies  peculiar  to  periodic  surfaces.  If,  for 
instance,  for  a  given  K,  rcn(K)  -  0  then  the  corresponding  P  wave  is  in  a  limbo:  it  is  neither  propagating  nor 
evanescent;  but  even  a  slight  alteration  of  the  frequency  could  make  it  of  either  type,  thereby  causing  a 
redistribution  of  energy  among  the  surviving  propagating  waves.  (Of  course,  this  redistribution  occurs  only  if 
any  of  tfe?:  refiracting  surfaces  is  periodically  corrugated.)  An  SV  wave  would  have  such  an  instability  if  its  on(K) 
«  0.  This  instability  is  called  the  Rayleigh-Wood  anomaly,  and  a  reasonable  description,  pertinent  to  the  present 
problem,  is  given  elsewhere  [Lakhtakia  et  al.,  1984b;  1985a].  Suffice  it  will  to  state  here  that  if  any  one 
member  of  the  four  sets  of  quantities,  rcn(K),  an(K),  rc'n(K),  and  a'n(K),  happens  to  vanish,  then  a  sharp 
Rayleigh-Wood  anomaly  can  be  observed  in  the  radiation  patterns.  As  a  numerical  example,  a  Rayleigh-Wood 
anomaly  occurs  in  the  P-wave  radiation  patterns  of  Fig.  5  at  9  «  ±9.3  deg  when  K+jfkpIsinSI)  -  ±kp,  which  is 
also  reflected  in  the  SV  radiation  patterns  at  9  —  ±5.4  deg.  Likewise,  at  -9  «  ±42.5  deg  when  K+2(kp|sin9|)  » 
±kp,  which  is  also  reflected  in  the  SV  radiation  patterns  at  -9  «  ±23.2  deg.  Also,  at  -9  -  ±27.8  deg, 
K^3(ks|sin9|)  «  ±kj,  showing  an  anomaly  in  the  SV  radiation  pattern.  There  are  no  such  anomalies  in  Fig.  7 


where  the  upper  surface  was  taken  to  be  planar.  Not  all  of  such  anomalies  are  very  prominent,  and  some  may 
have  indeed  been  missed  because  of  the  discretization  of  the  range  of  tc  for  ease  in  calculations. 


Finally,  the  last  identifiable  reason  for  some  of  the  various  anomalies  to  exist,  when  Tb  is  a  bimaterial 
interface,  is  the  excitation  of  a  Stoneley  wave  [Cagniard  et  al.,  1962;  Yamaguchi  and  Sato,  1955;  Phynney, 
1961]  on  Tb.  A  Stoneley  wave  is  a  non-dispersive  surface  wave  which  propagates  on  a  solid-solid  interface;  it 
decays  as  y‘^2  away  from  the  interface,  and  its  x-directed  wavenumber  can  have  two  roots  which  are  not 
necessarily  real  nor  on  the  proper  Riemann  surface.  If,  however,  for  a  K  :  |tc|  £  ji/L,  there  exists  a  real  k^k) 
which  satisfies  the  condition 

0  -  k„2  K2  "  (l/2)co2fp'-p)/(n'-n)]2  +  nn  an  [i^2  -  (l/2)ca2p’/(ti'-n)]2  + 

+  Jt'n  cr’n  [k„2  +  (l/2)co2p/((i'-|i)]2  ♦  [nn  an  +  n’n  a’J  (l/4)w2p  p’/(M'-l02  +  "n^'n^n  Kn2- 


then  a  Stoneley  wave  will  be  excited  for  that  k,  and  its  effect  can  be  observed  in  the  P-wave  radiation  pattern  at  9 
■  ±sin‘*(K/k'_ ),  and  in  the  SV-wave  radiation  pattern  at  9  -  ±sin‘*(ic/k’  ).  The  subject  of  Stoneley  waves  is,  and 
has  been,  in  itself  meritorious  of  extensive  attention,  for  which  reason  we  shall  not  discuss  it  any  further  here. 

Having  gone  through,  qualitatively,  the  various  reasons  for  the  occurrence  of  anomalies  in  the  radiation 
patterns,  the  characteristics  of  the  displacement  field  excited  on  the  traction-free  surface  Ta  remain  to  be 
investigated.  This  is  an  important  field,  very  often  being  measured,  particularly,  in  seismological  studies. 
Shown  in  Fig.  8  are  the  plots  of  the  two  components  of  ua  versus  ±x/L  when  the  line  source  is  located  at  Tp  * 
(0,0)  an''  the  normalized  frequency  kpL  =  7.5.  The  wave  velocities  in,  and  the  densities  of,  the  two  media  are  the 
same  as  in  Fig.  2;  a/L  -  d/L  -  1.0;  and  b/L  -  0.05.  In  Fig.  9,  these  same  quantities  are  re-plotted  but  for  b/L  « 
0.0  so  that  Ta  is  flat  The  highlight  nf  these  figures  is  the  extreme  variability  of  the  components  of  ug.  In 
particular,  when  Ta  is  flat,  then  bunching  of  data  seems  to  occur,  with  some  evidence  of  quasi-periodicity,  as 
may  be  seen  by  the  repetition  of  similar  (not  necessarily,  the  same)  features  in  Fig.  9.  In  Fig.  8,  on  the  other 
hand,  the  periodicity  of  Ta  creates  additional  anomalies  due  to  inter-modal  interactions;  hence,  the  data  variability 
is  enhanced,  and  the  bunching  is  diffused.  Incidently,  the  computations  of  ua  were  carried  out  by  a  63-point 
Gauss-Legendre  quadrature  integration  [Abramowitz  and  Stegun,  1972]  for  the  integral  over  k  :  |Kj  <>  rt/L  in  (19), 
and  considering  upto  n  «  ±5  for  each  K,  keeping,  of  course,  in  mind  that  the  computations  for  each  K  converged 
satisfactorily. 


are  present  all  spatial  frequencies  (tc  :  |tc|  £  <*>)  [Goodman,  1968],  having  the  same  amplitude  but  differing  in 
phase  when  starting  out  from  rp.  Therefore,  in  the  present  context,  there  is  no  likelihood  that  ua  will  eventually 
decay  down  uniformly  for  some  relatively  small  values  of  ±x/L  on  Ta,  as  can  be  observed  from  the  range  |x/L|  £ 
100.0  of  Figs.  8  and  9.  Parenthetically,  it  is  noted  here  that  were  the  line  source  to  be  an  uni-directional  radiator, 
radiating  out  both  P  and  SV  waves,  then  ua  would  decay  down  uniformly  to  some  mean  value  (zero  or 
otherwise)  at  some  reasonably  small  values  of  ±x/L.  This  conclusion  has  been  drawn  also  by  Fokkema  [1979] 
for  the  case  when  the  media  on  either  side  of  are  identical. 

The  last  data  set  computed  is  shown  in  Figs.  16  and  17  where  juaj|  and  lu^l.  respectively,  have  been 
plotted  against  the  normalized  frequency  kpL  for  some  selected  locations  ±x/L  on  Ta.  The  wave  velocities  in, 
and  the  densities  of,  the  two  media  involved  are  the  same  as  in  Fig.  2;  and  Ta  is  either  sinusoidally  corrugated 
(b/L  -  0.05)  or  flat  (b/L  -  0.0).  At  x  «  0,  uaj  is  identically  zero,  for  which  reason  the  latter  has  not  been  plotted 
in  Fig.  16.  From  these  graphs  it  is  observed  that  generally  ua  decreases  in  magnitude  as  |x/L|  and  kpL  become 
quite  large.  This  is  not  surprising,  and  for  the  following  reason:  at  any  location  r,  the  incident  Field  u*  of  Eq.  (9) 

possesses  variations  of  the  kind  |r  -  r^'1,  |r  -  rpf2.  |r  -  rp|'3 . When  r  is  far  away  from  the  source  location  r^ 

then  the  only  pertinent  variation  of  u1  is  |r  -  r^"1,  the  Field  magnitude  decaying  inversely  with  the  distance  |r  - 
rp|.  But  when  r  gets  closer  to  r^  then  the  higher  order  variations  also  become  significant,  and  the  Field  amplitude 
can  become  quite  large  [Lakhtakia  and  Iskander,  1983;  Lakhtakia  et  al.,  1984a].  This  means  that  although  the 
line  source  radiates  a  fixed  amount  of  power,  at  larger  distances  from  it,  its  radiated  field  can  be  quite  small, 
while  close  to  it,  the  Field  magnitude  can  be  very  high.  Furthermore,  the  near-zone  (i.e.,  the  zone  in  which  the 
higher-order  variations  of  the  incident  Field  are  significant)  has  a  larger  expanse  when  the  frequency  is  small,  and 
vice  versa.  Hence,  the  observation  just  made  from  Figs.  16  and  17. 

In  conclusion,  a  method  based  on  the  extended  boundary  condition  method  has  been  described  here  to 
compute  the  radiation  characteristics  of  an  elastodynamic,  isotropic  line  source  embedded  in  layered  media  with 
periodically  varying  interfaces.  Since  it  is  a  matrix  approach,  it  is  computationally  more  tractable  than  a 
procedure  involving  the  solution  of  integral  equations.  Although  the  speciFic  problem  solved  here  is  quite 
simple,  it  is  an  effective  example  of  a  wide  class  of  problems,  and  the  presented  approach  can  be  used  with  profit 
for  solving  them.  Because  of  the  specific  expansion  (17)  used  for  the  surface  field  on  Ta,  one  limitation  of  this 
method  is  that  the  maximum  slopes  of  the  corrugated  interfaces  may  not  exceed  0.448.  This  limitation  is  often 
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Fig.  2  Far-zone  P  wave  and  SV  wave  ( - )  radiated  amplitudes,  in  the  lower  half-space,  due  to  the 

P-wave  line  source  located  at  Tp  -  (0,0).  These  amplitudes  are  drawn  as  functions  of  8,  |0|  £  90  deg,  the  angle 
with  respect  to  the  -y  axis.  The  surface  Ta  is  specified  by  y  -  a  -  b  cos(2rcx/L)  while  is  the  plane  y  -  -  d. 
The  relevant  parameters  are  Cp  -  3500  m  s'*,  c$  -  6000  m  s'*,  p  -  2700  kg  m'3,  c’p  -  4500  ms*,  c'$  - 
8000  m  s'*  and  p’  -  3200  kg  m'3;  while  a/L  -  d/L  -  1.0  and  b/L  -  O.OS.The  normalized  frequency  (a)  kpL  - 
5.0,  (b)  kpL  -  7.5  and  (c)  l^L  -  10.0. 
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FI  t .  1 5  Same  as  Fig.  9,  eacept  Fq.  (31)  is  used  in  place  of  Eq.  ( 10)  for  the  incidenl  field. 
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Fig.  16  Computed  values  of  the  surface  displacement  component  uaj  as  function  of  the 


kpL  at  several  distances  x/L  along  the  sinusoidally  corrugated  Ta.  The  line  source  is  loc 


-  d/L  -  1.0;  b/L  -  0.05  ( - )  and  b/L  >  0.0  The  other  relevant  parameters  are 


6000  m  s'1,  p  -  2700  kg  m'3,  c’  -  4500  m  s'1,  c*  -  8000  m  s'1  and  p’  -  3200  kg  m' 


x/L  -  ±5.0;  (c)  x/L  -  ±10.0;  (d)  x/L  -  ±100.0. 
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SECTION  5 


The  Meckering  Earthquake  of  October  14,  1968. 
A  Possible  Downward  Propagating  Rupture. 

by 


Kristfn  S.  Vogfjord  and  Charles  A.  Langston 
Department  of  Geosciences 
The  Pennsylvania  State  University 
University  Park,  PA  16802 


ABSTRACT 


Average  source  parameters  of  the  1968  Meckering, 
Australia  earthquake  are  obtained  by  the  inversion  of  body 
waves.  The  objectives  of  the  inversion  are  the  elements  of 
the  moment  tensor  and  the  source  time  history.  An  optimum 
source  depth  of  3  km  is  determined,  but  because  of  source 
complexity  the  point  source  assumption  fails  and  the  moment 
tensor  obtained  at  that  depth  has  a  large  non-double-couple 
term,  CLVD  =  34  %.  The  source  parameters  of  the  major 
double-couple  are  :  strike  341°,  dip  37°,  rake  61°  and 
seismic  moment  8.2  x  1023  dyne-cm.  The  source  time  function 
is  of  approximately  4  seconds  duration,  with  a  long  rise 
time  and  a  sharp  fall-off.  The  fault  length  is  constrained 
on  the  surface  by  the  observed  surface  break  and  results 
from  vertical  displacement  modeling  suggest  a  width  of 
approximately  10  km  in  the  middle,  assuming  a  dip  of  37°. 
That  restricts  the  entire  faulted  area  to  lie  above  6  km 
depth.  Two  finite  fault  models  for  the  earthquake  are 
presented,  with  rupture  initiating  at  a  point  1)  near  the 
top  of  the  fault  and  2)  at  the  bottom  of  the  fault.  Both 
models  produce  similar  long  period  synthetics,  but  based  on 
the  short  period  waveforms,  model  1  is  favored.  It  is 
argued  that  such  a  rupture  process  is  the  most  reasonable  in 
this  cold  shield  region. 


1974  ;  Drummond,  1979).  Most  of  the  seismicity  in  the 
region  is  of  low  magnitude,  yet  two  of  Australia's  largest 
earthquakes  were  located  within  this  seismic  zone.  One 
occured  in  1941  and  was  located  in  the  northern  part  of  the 
zone,  the  other  was  the  Meckering  earthquake,  which  occured 
on  October  14,  1968. 

The  Meckering  earthquake  was  a  moderate-sized  intraplate 
event  with  a  surface  wave  magnitude  of  6.8  and  a  very 
shallow  source  depth,  1  kn,  as  determined  by  the  USCGS .  The 
epicenter  was  located  near  the  small  town  of  Meckering  which 
was  destroyed  in  the  earthquake  (Everingham  et  al.,  1969). 
This  event  is  of  special  interest  because  of  its  location  in 
Archaean  shield  within  the  South  West  Seismic  Zone.  Also  it 
is  one  of  a  few  intraplate  earthquakes  known  to  be 
accompanied  by  large  scale  surface  thrust  faulting.  The 
near  semicircular  appearance  of  the  fault  scarp  is  also 
quite  unusual.  The  surface  faulting  accompanying  the 
earthquake  and  its  aftershocks  extends  over  a  200  kma  area 
which  is  largely  covered  by  several  meters  of  sand  and  clay. 
Faulting  is  not  exposed  in  the  underlying  bedrock  which 
consists  mainly  of  deeply  weathered  granite  and  gneiss 
(Gordon,  1971).  The  largest  fault  is  the  Meckering  fault,  a 
37  km  long  arcuate  fault  scarp  (see  Figure  1)  concave  to  the 
east  and  dipping  about  40°  inward  (Gordon,  1971).  This 


fault  was  not  shown  on  any  structural  maps  of  the  region, 
prior  to  the  earthquake.  However,  Gordon  and  Lewis  (1980) 
report  findings  of  iron  stained  fault  breccia  on  the  fault, 
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and  during  the  mapping  of  the  fault  scarp  a  correlation  was 
noted  between  the  fault  trace  and  frequent  quartz  scatter, 
or  iron-rich  reddish  soil.  This  they  interpreted  as 
indications  of  a  preexisting  fault.  The  displacement  on  the 
fault  is  in  the  reverse-dextral  sense  with  the  eastern  side 
uplifted  and  thrust  over  the  western  side.  It  is  greatest 
near  the  middle  and  decreases  towards  the  ends.  The 
measured  throw  in  the  middle  is  about  2  m,  the  heave  about 
2.4  m,  and  the  dextral  movement  about  1.5  m  (Gordon,  1971). 

Stress  measurements  made  in  the  area  after  the 
earthquake,  at  depths  less  than  10  m  (Denham  et  al . ,  198°) 

indicate  a  high  level  of  horizontal  compressive  stress  in 
the  upper  crust,  with  the  average  direction  of  maximum 
principal  stress  about  N77E.  The  lowest  compress ional 
stress  was  measured  near  the  earthquake’s  epicenter, 
approximately  40  bars  and  increased  with  distance  away  from 
it,  up  to  approximately  230  bars  90  km  north  of  Meckering. 
This  large  horizontal  compressive  stress  is  the  probable 
cause  of  the  earthquake  and  assuming  that  the  difference  in 
stress  is  due  to  the  stress  relieved  by  the  earthquake,  it 
would  imply  a  stress  drop  of  about  100  bars.  The  dominance 
of  compression  observed  in  the  Australian  continent  can  be 
generated  by  plate  tectonic  forces  acting  on  the 
Indian-Australian  plate  (Cloetingh  and  Wortel,  1985). 
However,  the  compressional  direction  predicted  for  Western 
Australia  does  not  agree  with  the  observed  direction,  which 
may  therefore  be  caused  by  regional  effects.  The  proximity 


to  the  continental  margin,  which  has  experienced  several 
kilometers  of  normal  offset  along  the  Darling  Fault,  may  be 
an  influence. 

The  objective  of  this  study  is  to  gain  some  insight  into 
the  stress  condition  existing  in  the  upper  crust  of  the 
Precambrian  shield.  The  state  of  stress  can  be  deduced  from 
the  source  mechanism  of  the  Meckering  earthquake,  which  is 
determined  by  use  of  seismic  and  static  uplift  data.  A 
linearized  inversion  of  body  waves  is  performed  to  obtain  a 
moment  tensor  and  time  function  for  the  event,  then  the 
measured  ground  displacement  is  modeled  by  elastic 
dislocation  theory  to  obtain  some  information  about  fault 
dimensions,  dip  and  displacement.  Finally,  forward  modeling 
is  attempted  using  finite-fault  dislocation  models.  Results 
from  two  different  fault  models  are  presented. 

Previous  work  on  this  earthquake  includes  a  first  motion 
body  wave  study  by  Fitch  et  al.  (1973)  to  obtain  a 
double-couple  solution.  They  favored  reverse-sinistral 
motion  on  a  fault  plane  striking  332°  and  dipping  68°W,  with 
a  seismic  moment  of  6.1  x  10as  dyne-cm.  This  is  not 
consistent  with  the  dip  and  reverse-dextral  displacement 
indicated  by  the  surface  rupture.  A  small  foreshock  with  a 
different  mechanism  that  occured  about  3.5  seconds  before 
the  main  shock  may  have  affected  their  results  such  that 
first  arrivals  at  some  stations  were  incorrectly  chosen. 

This  study  offers  an  improvement  to  their  results  by 
providing  a  fecal  mechanism  that  is  in  better  agreement  with 


model  is  presented  in  Table  1.  Responses  from  73  P  rays  and 
18  SH  rays,  which  included  all  significant  arrivals  within 
the  first  30  seconds  were  calculated  and  summed  up.  The 
responses  were  then  convolved  with  the  WWSSN  15-100 
instrument  response,  and  an  attenuation  operator  with  t*  = 
1.0  for  P  waves  and  t*  =  4.0  for  S  waves.  The  source  time 
function  was  parameterized  as  a  series  of  weighted  boxcars, 
each  of  1  second  duration. 

Inversions  for  the  moment  tensor  and  six  time  function 
elements  were  performed  with  the  source  at  2  km  intervals  in 
the  depth  range  of  1-7  km.  Five  iterations  were  performed 
at  each  depth  but  convergence  was  usually  obtained  within 
three.  In  this  study,  the  parameters  were  not  weighted  and 
the  covariance  matrix  for  the  data  was  assumed  to  be  the 
identity  matrix.  A  cut-off  value  specifying  the  maximum 
allowable  variance  of  the  parameter  changes  was  set  to  be 
0.3.  With  that  choice,  all  parameters  were  well  resolved 
and  a  unique  solution  obtained.  The  condition  numbers  for 
the  inv<-  r3i  en  .  ware  between  60  and  100. 

The  final  moment  tensor  is  decomposed  into  a  major 
double-couple,  which  is  the  average  of  the  maximum  and 
minimum  principal  components  and  a  remainder,  which  is  the 
compensated  linear  vector  dipole,  or  CLVD .  The  relative 
size  of  the  CLVD  can  be  a  measure  of  the  complexity  of  the 
source,  although  a  large  CLVD  component  can  also  be  an 
indication  of  deficiencies  in  the  model.  Other  indicators 
of  the  performance  of  the  inversion  are  the  RMS  error,  which 


measures  the  fit  of  the  synthetics  to  the  data,  and  the 
least  squares  error  which  is  a  measure  of  the  quality  of  the 
inversion.  Figure  2  a)  shows  the  variation  of  the  CLVD  with 
depth.  It  is  quite  large  at  all  depths,  although  it 
decreases  rapidly  as  the  source  depth  approaches  the 
surface.  Figure  2  b)  shows  the  change  in  RMS  error  with 
source  depth  and  in  Figure  2  c>,  the  variation  of  the  least 
squares  error  with  depth  is  shown.  Although  resolution  is 
not  good,  especially  in  the  RMS  error,  they  both  have  minima 
at  3km  depth. 

Incorrect  assumptions  made  in  calculating  the  Green’s 
functions,  such  as  source  depth  and  the  value  of  t*  can  be 
compensated  by  changes  in  the  source  time  function,  to 
produce  the  same  synthetics.  An  error  in  t*  is  compensated 
by  the  length  of  the  source  function,  and  deviations  from 
the  correct  source  depth  add  complications  to  the  source 
function.  The  simplest  time  function  therefore  corresponds 
to  the  best  estimate  of  the  source  depth  (Christensen  and 
Ruff,  1985).  The  inversion  time  functions  for  various 
source  depths  are  plotted  in  Figure  3,  where  the  simplest 
time  functions  resjlt  from  the  3  km  and  5  km  source  depths, 
and  the  deepest  source  at  6.9  km  produces  a  periodically 
oscillating  time  function,  indicating  that  the  actual  source 
depth  has  been  exceeded.  These  results  combined  with  the 
minima  in  the  least  squares  and  the  RMS  error  indicate  a 
source  depth  of  approximately  3  km.  The  duration  of  the 
time  function  was  not  known  in  advance,  so  in  order  to  avoid 


orientations  thus  obtained  are  better  constrained  and  give  a 
near  vertical  tension  axis. 

The  large  non-double-couple  term  in  the  moment  tensor 
indicates  that  the  assumed  point  source  model  may  be 
inadequate  for  determining  the  source  mechanism,  and  a 
finite-sized  source  allowing  rupture  on  a  non-planar  fault 
might  have  to  be  incorporated  to  explain  the  radiation 
pattern  completely.  In  fact,  the  surface  rupture  strongly 
suggests  a  moment  tensor  with  a  changing  source  orientation. 
Lateral  heterogeneity  of  the  source  region  can  also  disturb 
the  radiation  pattern,  so  the  proximity  of  the  source  area 
to  the  shield  boundary  where  i«>ajor  lateral  discontinuities 
occur  is  also  of  some  concern.  Therefore,  the  assumption  of 
a  homogeneous  plane-layered  source  structure  is  likely  to  be 
the  cause  of  tome  error  in  the  solution,  especially  at 
stations  in  western  azimuths. 


Ill . 


surface  displacement 


The  surface  deformation  caused  by  the  earthquake  can  be 
used  to  place  some  constraints  on  the  fault’s  dimensions  and 
dip.  The  static  data  available  are  presented  in  Figure  6 
and  consist  of  throw  measurements  on  the  fault  scarp,  and 
uplift  measurements  obtained  from  relevelling  surveys  of 
pipelines  and  the  road  t^  Meckering  (Gordon,  1971).  Also 
shown  are  the  contours  of  equal  uplift,  which  for  lack  of 


resulting  moment  is  14  x  10as  dyne-cm  which  is  nearly  two 
times  the  value  obtained  from  the  body  wave  inversion.  This 
difference  suggests  that  roughly  half  of  the  energy  released 
in  the  earthquake  excited  the  body  waves  and  the  rest  may 
have  been  released  by  slow  slip  on  the  fault  after  the 
initial  rupture  phase,  and  or  by  the  aftershocks. 
Overestimated  fault  area  can  also  be  a  factor.  Using  a 
relation  from  Kanamori  and  Anderson  (1975)  for  the  stress 
drop  on  a  shallow  dip-slip  fault, 

Ac  *  m  U  /  W 

where  W  is  the  width  of  the  fault,  the  stress  drop  is  found 
to  be  :  Ao  =s  100  bars.  Fitch  et  al .  (  1973)  obtained  a 
similar  stress  drop,  although  their  estimate  of  displacement 
(1.5  m)  was  too  small  and  they  assumed  a  fault  width  of  only 
5  km,  which  is  half  the  width  obtained  in  this  study. 

IV.  Fii.  te  Fault  Models 

An  attempt  was  made  to  produce  matching  synthetics  with 
the  use  of  finite  fault  models.  The  method,  described  in 
Langston  (1978)  is  identical  to  the  method  used  to  calculate 
the  Green’s  functions  for  the  inversion,  except  for  the 
calculation  of  the  source  time  function.  A  propagating 
dislocation  source  with  a  given  rupture  velocity  on  a  finite 


fault  is  assumed.  The  dislocation  at  each  point  on  the 
fault  is  described  by  a  step  function,  delayed  in  time 
relative  to  the  distance  away  from  the  initation  of  rupture, 
and  the  resulting  source  time  function  is  the  integral  of 
the  dislocation  over  the  whole  fault  area.  The  source  time 
functions  thus  computed  differ  for  P  and  S  waves,  and  for 
up-  and  downgoing  rays  from  the  initial  source.  Depending 
on  the  directivity,  some  variation  with  distance  and  azimuth 
may  also  occur  in  the  source  functions.  In  addition  to  the 
long  period  P  waveforms  used  in  the  inversion,  short  period 
P  waves  were  also  included  in  this  part  of  the  study. 

Because  of  the  foreshock  3.5  seconds  before  the  main 
shock  and  the  shallow  source  depth,  extreme  interference 
makes  it  impossible  to  identify  distinct  phases  in  the  short 
period  seismograms;  the  P  phase  of  the  main  shock  is  buried 
in  the  coda  of  the  foreshock  and  the  P,  pP  and  sP  phases  all 
arrive  within  a  short  time  interval,  thus  causing  the 
interference  and  making  amplitudes  small  at  the  beginning  of 
the  records.  However,  by  careful  alignment  in  time  of  the 
long  period  and  the  short  period  waveforms,  the  arrival  of 
the  main  shock  in  the  short  period  record  can  be  inferred 
from  its  arrival  in  the  long  period  record.  This  is 
demonstrated  in  Figure  9,  where  the  first  arrow  indicates 
the  arrival  of  the  foreshock  and  the  second  arrow  indicates 
the  main  shock  arrival. 

Starting  out  with  the  fault  parameters  of  the  major 
double-couple  from  the  inversion  and  a  fault  width  of  9.5  km 


in  the  middle,  as  determined  by  the  surface  uplift,  best 
fitting  synthetics  were  calculated  by  trial  and  error 
procedure.  Synthetics  were  made  assuming  either  a  rupture 
initiation  at  a  point,  or  a  propagating  finite  line  source. 
The  parameters  allowed  to  vary  were  the  fault  dimensions 
(except  for  its  width  in  the  middle),  the  rupture  velocity, 
and  the  initial  source  depth.  All  these  parameters  can 
trade  off  to  produce  the  same  synthetics  for  various 
different  choices  of  parameters.  The  best  fitting  long 
period  and  short  period  synthetics  were  obtained  with 
rupture  initating  at  a  point  at  1  km  depth  and  propagating 
radially  with  a  rupture  velocity  of  2.0  km/sec.  The 
orientation  of  the  near  semicircular  fault  plane,  which 
extends  down  to  5.7  km  depth  is  shown  in  Figure  10.  The 
long  period  and  short  period  waveforms  produced  by  this 
model  are  plotted  in  Figures  11  and  12,  respectively.  The 
overall  shape  of  the  long  period  waveforms  is  fairly  well 
fit  at  all  stations,  although  the  first  motions  at.  the 
stations  in  western  azimuths  are  not  perfectly  matched.  The 
short  period  waveforms  are  also  surprisingly  well  fit, 
considering  that  such  a  simplified  model  of  the  fault  is 
used . 

The  dominant  factors  in  shaping  the  short  period 
synthetics  are  the  source  time  functions  of  the  up-  and 
downgoing  rays.  Due  to  the  slowly  increasing  fault  area  and 
then  the  abrupt  stopping  of  rupture,  the  source  time 
functions  have  a  long  rise  time  and  a  sharp  fall-off.  That 


causes  the  receiver  responses  to  have  the  largest  amplitudes 
at  the  end.  When  the  responses  due  to  the  P,  pP  and  sP 
phases  are  added  up  they  produce  the  large  negative 
amplitudes  denoted  by  A  and  B  in  Figure  12,  where  A  is  the 
backswing  of  the  P  response,  and  B  is  due  to  the  end  of  the 
pP  or  sP  response,  depending  on  which  is  the  dominant  one. 

In  eastern  azimuths,  at  stations  GUA,  AFI ,  and  RAR,  the  sP 
phase  is  the  largest,  but  at  all  other  stations  pP  has  the 
largest  amplitude .  The  arrivals  of  these  phases  are  also 
denoted  in  the  figure.  The  main  features  in  the  short 
period  synthetics  that  simulate  the  seismograms  so  well  are 
therefore  caused  by  the  sudden  stopping  of  rupture,  and  for 
such  a  shallow  source  depth  the  P,  pP  and  sP  phases  all 
arrive  within  the  first  second,  causing  destructive 
interference  and  very  small  amplitudes  in  the  first  few 
seconds.  The  observed  short  period  seismograms  also  have 
small  amplitudes  within  the  first  seconds,  but  they  are 
contaminated  by  the  coda  of  the  foreshock,  making  it 
impossible  to  identify  specific  phases.  A  few  of  the  short 
period  waveforms  are  not  well  matched  by  the  synthetics. 

These  are  generally  richer  in  high  frequencies,  which  can  be 
explained  by  a  lesser  attenuation  along  these  raypaths. 

Small  changes  in  the  orientation  and  dip  of  the  fault 
might  improve  the  first  motion  fits  in  the  long  period 
waveforms  and  would  not  cause  any  major  changes  in  the  short 
period  waveforms,  as  long  as  the  source  depth  and  fault 
dimensions  are  kept  constant. 


The  fault  area  of  the  model  is  144  km  ,  and  the  seismic 
moment  is  8.2  x  102S  dyne-cm.  That  gives  a  displacement  of 
1.6  m  on  the  fault,  which  is  less  than  the  displacement 
inferred  from  the  static  data. 

A  second  model  was  also  tried  to  see  if  matching 
synthetics  could  be  made  with  rupture  initiating  at  depth 
and  propagating  upwards.  The  fault  was  made  to  approximate 
the  surface  break  more  closely  and  the  rupture  was  started 
at  a  point  at  5.2  km  depth,  at  the  bottom  of  the  fault,  and 
propagated  radially.  The  orientation  of  this  fault  model  is 
shown  in  Figure  13.  The  fault  was  broken  into  three 
sections  with  equal  dips  of  30°  but  different  strikes.  The 
direction  of  slip  was  the  same  on  all  sections.  By  trial 
and  error,  best  fits  were  obtained  with  a  rupture  velocity 
of  2.7  km/sec  and  a  30  %  higher  moment  on  the  middle  section 
tha.n  on  the  other  two.  The  total  moment  was  8.0  x  1023 
dyne-cm  and  the  total  area  was  194  km2.  The  corresponding 
slip  on  each  section  is  1.8  m  in  the  middle  and  1.0  m  on  the 
other  two.  The  long  period  waveforms  produced  by  this  model 
are  plotted  in  Figure  14.  They  fit  the  seismograms  fairly 
well  a.  ?  the  fit  is  as  good  as  for  the  previous  model. 
However,  the  short  period  waveforms,  which  are  plotted  in 
Figure  15  do  not  fit  the  seismograms  very  well,  the  main 
reason  being  that  the  time  functions  for  the  down  going  rays 
are  too  long,  so  that  the  backswing  from  the  P  arrival  is 
annihilated  by  the  backswing  from  the  surface  reflected 
phases,  pP  and  sP.  The  time  intervals  between  the  arrivals 


zone,  in  the  ductile  region,  the  whole  rock  matrix  yields  by 
aseismic  creep  and  the  shear  strength  decreases  drastically 
with  depth. 

Compiled  stress  measurements  from  various  regions 
( McGarr ,  1980;  Bamford,  1976)  indicate  a  linear  increase  in 
shear  stress  with  depth,  but  the  gradient  is  less  than  that 
of  the  rock’s  static  frictional  strength.  Accordingly,  the 
shear  stress  is  generally  less  than  the  rock  strength  over 
all  of  the  brittle  region.  In  the  ductile  region  however, 
shear  stress  and  rock  strength  will  eventually  become  equal 
and  aseismic  slip  can  occur.  In  areas  of  medium  to  high 
heat  flow  this  occurs  at  shallow  depths  and  the  ductile 
deformation  can  build  up  enough  stress  in  the  transition 
zone  above  to  allow  rupture  to  initiate  on  a  crack  and 
continue  upwards  to  cause  an  earthquake.  The  rupture  stops 
vhent  the  shear  stress  becomes  smaller  than  the  dynamic 
frictional  strength  of  the  rock.  This  could  be  the  case  for 
the  aforementioned  events,  which  all  initiated  in  or  above 
the  transition  zone  and  ruptured  upwards.  On  the  other 
hand,  in  a  cold  region  such  as  South  West  Australia,  by 
assuming  a  dry  lower  crust  (Meissner  and  Strehlau,  1982)  the 
transition  zone  is  at  approximately  30  km  depth.  No 
earthquakes  have  been  observed  at  that  depth  in  the  region, 
implying  a  shear  stress  in  the  lower  crust  that  is  too  low 
to  overcome  the  static  frictional  strength  and  cause 
rupture.  However,  the  measured  shear  stress  at  the  surface 
is  quite  high,  about  100  bars  (Denham  et  al . ,  1980)  and 
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assuming  the  same  linear  gradient  as  before,  the  shear 
stress  can  exceed  the  frictional  strength  in  a  limited  depth 
interval  just  below  the  surface.  There  the  stress 
difference  should  also  be  greatest,  and  decrease  with  depth. 
Accordingly,  rupture  is  most  likely  to  initiate  near  the 
surface  and  propagate  to  some  depth  before  terminating.  This 
model  therefore  indicates  that  in  the  cold  shield  of  Western 
Australia,  the  brittle-ductile  transition  zone  is  far  below 
the  fault  area,  which  lies  entirely  in  the  brittle  region, 
and  only  at  shallow  depths  is  the  shear  stress  sufficiently 
high  to  cause  rupture. 

The  horizontal  compression  of  the  crust  in  Western 
Australia  does  not  explain  the  existence  or  the  location  of 
the  Seismic  Zone.  The  driving  force  for  the  seismicity 
could  be  isostatic  movements  of  the  Yilgarn  Block,  perhaps 
due  to  the  thickening  of  the  dense  basal  layer  towards  the 
shield  boundary.  This  vertical  movement  could  cause 
compressive  stresses  over  a  limited  area  at  the  surface  that 
would  add  just  enough  to  the  regional  compressive  stress  to 
trigger  earthquakes. 


VI .  Conclusions 

A  moment  tensor  solution  that  describes  the  average 
fault  mechanism  of  the  Meckering  earthquake,  has  been 
obtained.  The  horizontal,  east-west  direction  of  the 
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Table  1.  -  Crustal  Model 

VP  ( km/sec )  V,(km/sec)  Density  (g/cmM  Thickness  (km) 


6.13 

3.54 

2.78 

7.0 

6.70 

3.87 

2.94 

7.0 
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Figure  5.  Observed  (above)  and  synthetic  (below)  SH 

waveforms  for  a  3  km  source  depth.  Waveforms 
are  normalized  to  their  maximum  amplitude,  which 
is  indicated  to  the  right  of  each  trace.  The 
inversion  time  windows  are  indicated  by  the 
arrows.  Also  shown  are  lower  hemisphere  equal - 
area  projections  of  the  moment  tensor  nodal 
surfaces  (solid  lines),  the  nodal  planes  of  its 
major  double-couple  (dashed  lines) ,  and  the 
location  and  polarities  of  SH  wave  first  motions 
(*  "clockwise  ;  0  "nodal  ;  -  "counterclockwise) 
The  pressure  (P),  tension  (T)  and  intermediate 
(I)  axes  ere  also  indicated. 
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Figure  6.  Map  showing  vertical  displacement  on  the  fault 
scarp  and  uplift  based  on  relevelling  of  a 
pipeline  and  the  road  from  Northam  to  Meckerlng 
(uplift  in  meters) .  The  uplift  contours  are 
estimated  away  from  the  scarp  and  the  relevelled 
lines  (Modified  from  Gordon,  1971). 
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DISTANCE  FROM  FAULT  SCARP  (km) 


Vertical  displacement  of  a  47*  dipping  fault 
(x'b),  and  a  27*  dipping  fault  (circles). 
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Summary 


Source  parameters  for  five  moderate  African  intraplate 
events  were  obtained  from  body  wave  inversion  for  the  moment 
tensor.  Parameters  determined  for  the  events  are  as  follows: 


location 

date 

Mb 

strike 

dip 

rake 

depth 

Ethiopia 

4/05/69 

5.8 

327.1 

69.2 

-45.0 

4 

Congo 

3/20/6U 

F  3 

44.6 

44.9 

-81.5 

10 

Zambia 

12/2/68 

5.9 

20.3 

40.2 

264.8 

6 

Zambia 

5/15/68 

5.7 

49.1 

39.8 

263.2 

28 

S.  Africa 

9/29/69 

5.6 

123.5 

89.4 

-0.58 

4 

The  Ethiopia 

and  Congo 

events 

are  both 

normal 

faults  as 

would 

be  expected  for  their  rift  zone  locations.  The  5/15/68  Zambia 
event  is  unusual  since  the  inversion  for  this  event  yielded  a 
high  angle  normal  fault  with  a  well  constrained  focal  depth, 
from  P-pP  times,  of  28km.  The  depth  and  focal  mechanism  of 
this  event  along  with  the  mechanism  for  the  12/2/68  Zambia 
event  suggest  the  possibility  of  southwest  rift  extension. 
Additionally,  the  May  15th  event  may  be  located  well  below 
the  brittle/ductile  transition  zone  which  implies  that  strain 
rates  in  the  area  are  higher  than  normal.  Teleseismic 
vertical  P  waves  from  the  South  Africa  event  exhibit 
distortions  that  may  be  attributed  to  lateral  inhomogeneities 
in  the  form  of  a  coastward  thinning  crust  at  the  source. 
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Introduction 

The  East  African  rift  system  is  one  of  the  most 
intriguing  tectonic  features  in  the  world.  It  was  one  of  the 
first  tectonic  features  attributed  to  the  extension  of  the 
earth’s  crust  (M=Kenzie  et  al.  1970).  Minster  et  al. 
(1974,1973)  have  not,  however,  been  able  to  reconcile 
extension  of  the  East  African  rift  with  kinematic  plate 
models  in  their  global  plate  modeling  studies.  The  stress 
regimes  associated  with  this  rifting  have  been  studied  by 
many  authors.  McKenzie  et  al.  (1970),  Maasha  &  Molnar 
(1972),  Fairhead  &  Girdler.  (1971),  and  this  paper,  attempt  to 
determine  stress  regimes  and  subsequent  plate  motions  by 
studying  the  focal  mechanisms  of  large  magnitude  (  >  5.5) 
earthquakes  of  the  Eastern  Rift  system.  Previous  authors 
used  P  wave  first  motion  from  long  and  short-period  P  records 
to  determine  source  mechanisms.  This  study,  conversely,  uses 
body  wave  inversion  for  the  moment  tensor,  and  extracts 
source  parameters  from  the  moment  tensor  elements  assuming  a 
deviatoric  point  source.  Inversion  is  advantageous  in  that 
it  allows  for  more  accurate  and  reliable  mechanisms  in 
addition  to  yielding  more  information  about  the  source.  The 
inversion  matches  the  amplitude  of  each  phase  and  allows  for 
the  determination  of  a  source  time  function  and  focal  depth. 
Substantial  information  about  source  crustal  structure  can 
also  be  obtained  by  matching  source  area  crustal  phases  in 
the  observed  data. 


The  East  African  rift  system,  at  a  glance,  has  two 
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Data  Analysis 

Long-period  WWSSN  P  and  S  wave  data  were  used  in  this 
study.  Vertical  P,  NS  and  EW  S  wave  seismograms  were  photo 
copied  from  70  mm  film  chips,  digitized  along  the  top  and 
bottom,  and  averaged  to  insure  smooth  wave  forms.  The 
digitized  data  were  linearly  interpolated  with  a  quarter 
second  sampling  interval  to  yield  the  digital  data,  The  NS 
and  EW  data  were  further  processed  by  rotating  them  into 
their  theoretical  back  azimuth  to  obtain  the  radial  and 
tangential  components  of  displacement. 

The  Green’ s  functions  needed  for  the  inversion  were 
computed  using  shallow  dislocation  source  theory  of  Langston 
&  Helmberger  (1975).  For  a  detailed  discussion  on  inversion 
the  reader  is  referred  to  papers  by  Wiggins  (1972),  Jackson 
(1972),  Langston  (1981),  and  Barker  &  Langston  (1981).  Figure 
2  compares  P  wave  lower  hemisphere  focal  spheres  obtained  by 
inversion  to  those  obtained  by  previous  authors.  Table  1 
gives  our  moment  tensor  elements,  source  depths,  source  time 
functions,  and  error  estimates. 

Inversion  error  was  determined  by  looking  at  the 
variance  of  the  parameter  changes  ]  and  the  resolution 

matrix.  The  (tables  1.1  &  1.2)  are  highly  dependent  on 

the  eigenvalues  with  small  eigenvalues  yielding  high  L&p  and 
a  highly  singular  inverse  opperator.  Our  were  computed 
assuming  *0%  error  in  the  data.  Because  of  the  importance  of 
the  eigenvalues,  we  have  included  the  condition  number, 
quotient  of  the  highest  and  lowest  eigenvalues,  in  table  1.3. 
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All  our  resolution  matricies  were  identity  matricies. 

Figure  3  illustrates  how  the  least  square  fit  between 
the  data  and  synthetics  varied  with  source  depth  for  each 
event;  our  final  focal  depths  were  those  corresponding  to  the 
minimum  least  square  error  each  event.  Schematics  of  the 
source  time  function  are  included  with  the  diagram  for  the 
Congo  event  to  show  how  both  the  time  function  simplicity  and 
the  least  square  fit  were  used  to  determine  the  best  focal 
depth  for  this  event.  These  time  functions  also  illustrate 
the  trade  off  between  source  depth  and  source  time  function 
duration  (each  bar  represents  a  one  second  boxcar),  with  the 
deeper  sources  having  shorter  time  functions. 

Ethiopia  -  April  5.  1969 

McKenzie  et  al.(1970)  were  the  first  to  determine  a 
mechanism  for  this  event.  Their  solution  has  fault  planes 
that  strike  333« ,  dip  90°,  and  strike  63«  ,  dip  70°SE.  They 
concluded  that  the  event  lies  on  a  dextral  transform  fault 
that  joins  the  East  African  rift  valley  and  the  Gulf  of 
Tad jura  because  of  observed  ground  displacement.  Fairhead  & 
Girdler  (1971),  conversely,  choose  the  63°  striking  plane 
because  of  reported  ENE  trending  af terschocks ;  this  choice  is 
consistent  with  the  rotation  of  the  Danakil  horst  away  from 
the  Nubian  Plate.  These  two  interpretations  illustrate  the 
opposing  views  of  the  tectonic  situation  in  the  Afar  region. 

Our  fault  planes  are  slightly  different  then  those 
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obtained  by  McKenzie  et  al.  with  strike  327°,  dip  69°NE,  and 
strike  77°,  dip  49°SSE  (figures  4a  &  b).  The  strikes  of  our 
planes  are  close  to  those  of  McKenzie  et  al.,  but  the  actual 
motion  along  these  planes  is  different.  We  arrived  at  a 
focal  depth  of  4  km  using  a  modified  —  thinned  by  9  km  — 
Gumper  &  Pomeroy  (1970)  crustal  model. 


Republic  of  Congo  -  March  20,  1966 

This  event  has  been  the  focus  of  both  first  motion  and 
aftershock  studies.  Sykes  (1967)  used  the  first  motion  of  P 
and  PKP  phases  from  long-period  records,  or  short-period 
where  no  long-period  data  were  available,  to  determine  a 
mechanism.  His  fault  planes  strike  14°,  dip  40°  W,  and 
strike  32°,  dip  52°  E.  Wohlenberg  (1966)  reported  a  fault 
with  a  40  km  extent  and  30-40  cm  of  vertical  displacement 
with  the  western  side  moving  up  relative  to  the  eastern  ^ide 
—  no  horizontal  displacement  was  observed.  Loupenkine  (1966) 
reported  a  15-20  km  rupture  striking  NNE  with  a  throw  of  2 
meters  --  the  down  throw  being  to  the  east.  Based  on  these 
reports,  Sykes  selected  the  plane  striking  14°  as  the  primary 
node.  His  axis  of  maximum  tension,  subsequently,  lies  nearly 
perpendicular  to  the  strike  of  the  rift.  No  attempt  was  made 
to  obtain  a  depth  of  focus. 

First  motion  P  data  alone  is  insufficient  to  constrain 
the  fault  planes  for  this  event  —  this  is  also  true  for  both 
Zambia  events.  Observational  data  (ie:  ground  displacement) 
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by  using  only  the  P  wave  data  illustrates  some  of  the 
problems  that  can  arise  with  inversion;  it  also  demonstrates 
the  trade  off  between  source  depth  and  time  function 
duration. 


Zambia  -  Mav  15,  1968 

Fairhead  &  Girdler  (1971)  obtained  a  fault  plane 
solution  using  first  motion  of  long-period  P  and  PKP  phases; 
they  obtained  a  normal  fault  striking  NNE.  No  surface  effect 
studies  are  available,  but  this  solution  is  in  good  agreement 
with  north-east  trending  faults  in  the  area.  But  because  the 
fault  did  not  rupture  the  surface,  and  because  the  P  wave 
first  motions  put  almost  no  constraints  on  the  fault  planes, 
Fairhead  &  Girdler’ s  mechanism  is  somewhat  questionable.  No 
attempt  was  made  to  determine  a  depth  of  focus.  Their 
mechanism,  nonetheless,  substantiates  a  WNW-ESE  extensional 
field  for  this  region. 

This  event  may  well  be  the  most  significant  of  the 
study.  Inversion  yielded  a  nearly  pure  dip  slip  fault  with  a 
28  km  depth  of  focus.  Teleseismic  P  waves  from  this  event 
show  that  it  occurred  at  depth;  similarly,  the  short-period 
records  are  a  textbook  example  of  how  pP-P  intervals  can  be 
used  to  determine  focal  depths  (figure  6).  While  this  event 
is  shallow  in  seismological  terms,  its  depth  is  3-7  times 
that  of  the  other  4  events  in  this  study.  Inversion  yielded 
a  high  angle,  nearly  pure  dip  slip  normal  fault  striking  at 
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either  49°  (dip  40°SE)  or  238°  (dip  51°NW),  with  a  4  second 
source  time  function  (figures  7a  &  h).  The  full  36  km  Gumper 
&  Pomeroy  crustal  model  was  used  in  the  inversion  and 
synthetic  computations  for  this  event. 


Zambia  -  December  2,  1968 

Fairhead  &  Girdler  (1971)  used  first  motions  of  P  and 
PKP  phases  to  obtain  a  mechanism  for  this  event  also.  Their 
two  planes  strike  27® ,  dip  38°ESE,  and  strike  57° ,  dip  55°NW. 
Their  normal  fault  mechanism  has  a  strike  similar  to  the  May 
15,  1968  Zambia  event  which  occurred  approximately  350  km  to 
the  south-east.  There  were  no  detectable  aftershocks  or 
observable  ground  displacement,  and  no  attempt  was  made  to 
obtain  a  depth  of  focus.  As  with  the  previous  event,  this 
mechanism  lacks  credibility  because  of  poor  first  motion 
constraints  and  a  lack  of  surface  rupture. 

Inversion  yielded  a  nearly  pure  normal  fault  (figures  8a 
&  b)  with  planes  striking  rr > ,  (dip  40°ESE)  or  207®,  (dip 
50°WNW) .  A  focal  depth  of  6  km  and  a  4  second  time  function 
yielded  the  lowest  least  square  error.  The  strike  of  this 
event  is  similar  to  that  of  the  May  15,  1968  Zambia  event; 
both  Zambia  events,  therefore,  provide  support  for  a  NWN-ESE 
extensional  field  for  south  central  Africa.  The  full  36  km 
thick  Gumper  fie  Pomeroy  crustal  model  was  used  in  inversion 
and  synthetic  computations  for  this  event  also. 
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accurate  solution.  This  was,  however,  not  the  case. 
Inversions  run  using  solely  P  data  resulted  in  mechanisms 
that  did  not  agree  with  P  first  motion  data.  Inversions  using 
only  SH  data,  conversely,  yielded  solutions  that  were 
consistent  with  both  P  and  SH  first  motions.  The  final 
solution  was  obtained  using  both  P  and  SH  waveforms  and  is 
consistent  with  all  first  motion  data.  This  problem, 
nonetheless,  was  very  puzzling  and  may  be  the  result  of 
complexities  in  the  crustal  structure  near  the  source. 

Langston  (1977)  notes  that  teleseismic  P  waves  radiated 
from  sources  whose  radiation  pattern  varies  rapidly  with 
azimuth  (ie:  strike  slip  faults)  can  be  greatly  distorted  by 
dipping  structure  at  the  source  —  this  is  particularly  true 
for  stations  located  near  the  nodal  surfaces.  We  assumed 
this  to  be  the  cause  of  poor  P  wave  inversion  results  and 
poor  P  wave  fits.  Synthetics  were  computed  using  Langston’s 
dipping  interface  algorithm  to  explore  the  effects  on  the 
synthetic  waveforms;  results  will  be  discussed  in  a  later 
section . 

Shudofsky  (1985)  used  body  wave  modeling  of  P  waves  in  a 
45°  sector  [342° -27° ]  to  determine  a  30  km  focal  depth  for 
this  event,  and  notes  similar  problems  in  his  waveform 
modeling.  He  concurs  with  the  nearly  pure  strike-slip 
mechanisms  obtained  in  other  studies,  but  comments  that  he 
could  not  model  the  complex  body  waves  using  the  simple 
source  properties  assumed  for  other  African  events.  He 
attributes  the  complex  wave  forms  to  source  directivity, 
multiple  sources,  or  unusual  source  structure. 
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Picking  first  motions  for  this  avent  was  difficult.  Some 
of  the  long-period  P  records  were  low  amplitude  and  had  very 
poor  signal  to  noise  ratios  which  mads  first  motion  picks 
unreliable.  Short-period  records,  when  they  were  available, 
were  used  to  determine  the  arrival  times  for  both  the  P  and  S 
waveforms.  Similarly,  long-period  records  from  stations  KOD, 
NDI,  and  SNG  (figure  9),  in  particular,  have  "arrivals"  that 
may  be  interpreted  as  P  and  sP.  Their  separation  suggest  a  9 
km  source  depth  which  is  not  substantiated  by  the  long-period 
SH  data  nor  the  short-period  P  or  S  data.  A  9  km  source  is, 
therefore,  ruled  out  in  favor  of  the  shallower  focus 
determined  from  the  short-period  records  and  the  least  square 
error  from  the  inversion. 

PiscvtgsiQn 


The  two  major  tectonic  questions  we  will  address  are: 
what  is  the  state  of  stress  at  various  locations,  and  how  far 
does  the  rift  system  extend  through  Africa.  We  have  divided 
our  study  area  into  3  sections  (figure  1):  section  #1  - 
northern  Afar  (figure  11),  sections  #2a  &  #2b  -  central 
Africa  (figures  13a  &  13b  respectively),  and  section  #3  - 
South  Africa  (figure  14). 


Afar  Region 


Two  plate  motion  configurations  have  been  postulated  for 
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moving  NNW  trending  normal  faults  are  a  better  indicator  of 
which  node  is  the  primary,  and  which  the  auxiliary.  Our 
preferred  mechanism  for  the  April  5,  1969  Ethiopia  event 
strikes  327®  and  dips  69°NE  (rake  -45® ) . 

Several  crustal  models  were  tested  for  validity  during 
the  inversion  for  each  event.  Griffiths  et  al .  (1971) 
proposed  a  20  km  two  layer  model  similar  to  models  for 
Iceland.  Bonjer  et  al.  (1971)  also  have  a  two  layer  model, 
with  the  crust  extending  to  a  depth  of  40  km.  Gumper  and 
Pomeroy  (1970)  suggested  a  3  layer,  36  km  thick  crust. 

All  three  models  and  two  additional  models, 
modifications  of  the  Gumper  &  Pomeroy  model,  were  used  in 
synthetic  waveform  computations  (figure  12).  A  thinned,  by  9 
km,  Gumper  &  Pomeroy  model  was  required  to  match  the 
[crustal]  reverberation  timing  of  the  observed  data  to  that 
of  the  synthetics,  particularly  for  waveforms  of  the  type 
observed  at  NDI  and  CHG.  The  27  km  thick  [model  a]  crust 
used  for  this  event  is  consistent  with  plate  tectonic  crustal 
thinning  assumptions  for  extensional  regimes. 

The  Afar  region  is  comprised  of  mantle  material  ejected 
during  the  separation  of  Arabia  and  Africa;  aeromagnetics 
show  that  this  accretion  occurred  at  a  very  slow  rate 
(McKenzie  et  al.  1970).  The  Ethiopia  event  demonstrates 
that  NNE-SSW  extensional  forces  still  exist  in  the  region  and 
subsequently,  that  the  Afar  region  is  still  growing  with  the 
Red  Sea.  The  low  level  of  both  seismic  and  volcanic  activity 
in  the  area  suggests  that  the  process  may  have  slowed,  but 
there  is  no  evidence  to  suggest  that  it  has  completely 
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stopped . 


Central  Africa 

The  second  area  we  examined  includes  regions  2a  and  2b 
of  figure  1.  Area  2a  (figure  13a)  encloses  the  Congo  event 
of  March  20,  1966;  2b  (figure  13b)  encloses  both  Zambia 
events  --  May  15,  1968  and  December  2,  1968.  The  Congo  event 
is  located  within  the  western  branch  of  the  rift  system, 
while  the  Zambia  events  are  located  in  an  area  that  may 
demonstrate  southward  continuation  of  the  rift.  The  Congo 
event,  which  occurred  in  the  Ruwenzori  massif  between  Lakes 
Albert  and  Edward  (figure  13a),  will  be  discussed  first. 

Two  independent  micro  earthquake  studies  (Rykounov  et 
al.  1972;  Maasha  1975)  have  shown  the  Ruwenzori  region  to 
be  one  of  the  most  seismically  active  areas  in  the  rift 
system.  This  area  of  the  rift  is  often  used  to  show 
similarities  between  the  mid-oceanic  spreading  centers  and 
the  Eastern  rift  system  (Baker  et  al.  1972).  Volcanism, 
seismic  velocities,  attenuation,  gravity  anomalies,  and 
seismicity  suggest  or  are  consistent  with  the  presence  of 
lateral  discontinuities  in  the  lithospheric  mantle  (Maasha  & 
Molnar  1972).  Micro-earthquake  data  show  a  prevalence  of  E-W 
tensional  stresses  (Maasha  1975).  This  information  shows 
that  the  March  20,  1966  Congo  event  occurred  in  a  region  that 
is  under  E-W  extensional  stress. 

Observational  data  can  be  used  to  determine  which  plane 
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is  the  primary  node.  While  the  two  reports  of  observable 
fault  displacement  (Loupenkine  1966;  Wohlenberg  1966) 
disagree  as  to  the  actual  extent  of  faulting,  both  agree  that 
the  western  side  of  the  fault  moved  up  relative  to  the 
eastern  side.  We,  therefore,  have  chosen  the  node  which 
strikes  45°  and  dips  45° SE  as  the  primary  node.  The  nearly 
pure  dip  slip  nature  of  this  fault  combined  with  its  10  km 
depth  show  that  the  crust  in  this  area  is  subject  to 
extensional  forces. 

The  26  km  thinned  Gumper  &  Pomeroy  crustal  model  was 
used  for  this  event.  As  with  the  Ethiopia  event,  it  is  not 
unreasonable  to  have  an  abnormally  thin  crust  in  this  area 
because  of  suspected  spreading  and  thinning.  In  any  case, 
crustal  thickness  has  little  effect  on  these  waveforms. 

The  two  Zambia  events  are  located  in  a  region  (figure 
13b)  that  is  under  scrutiny  for  some  proof  of  southward  rift 
extension.  Several  authors  (Reeves  &  Hutchins  1975;  Scholz 
et  al.  1976;  Chapman  &  Pollack  1977)  have  expressed  their 
views  on  the  possibility  of  rift  growth  that  was  first 
proposed  by  Fairhead  &  Girdler  (1969)  who  postulated  that  the 
western  branch  of  the  rift  followed  Lake  Tanganyika  and  then 
continued  south  through  Lake  Mweru.  They  suggested  that  rift 
extension  could  continue  through  Zambia  and  then  as  far  south 
as  24° S,  roughly  along  26° E. 

There  is  a  growing  body  of  evidence  supporting  the  rift 
extension  theory.  To  the  south  of  the  Zambia  events,  in 
Botswana,  Reeves  &  Hutchins  (1975)  detected  a  linear  feature 
they  named  the  Kalahari  Seismicity  Axis  [KSA] .  The  location 
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of  this  feature  was  determined  by  epicenter  data  and 
supported  by  photogeological  features.  The  KSA  strikes  NE-SW 
and  divides  the  Bouguer  anomaly  map  of  the  area  into  areas  of 
distinctly  different  textures  --  a  Bouguer  anomaly  gradient 
follows  much  of  the  axis.  The  authors  used  this  zone  to 
demonstrate  the  existence  of  tectonic  instability  in  the 
heart  of  the  continental  plate. 

Scholz  et  al.  (1976)  note  a  zone  of  seismicity  that  runs 
south-west  from  Lake  Tanganyika  through  Lake  Mweru  and  then 
south  through  central  Zambia  to  northern  Botswana.  Micro 
seismicity  in  this  area  is  higher  then  that  in  Ethiopia,  and 
comparable  to  that  observed  along  the  main  rift  valleys  of 
East  Africa.  Active  normal  faulting  on  north-easterly 
striking  faults  provides  evidence  of  a  very  recent  onset  of 
rifting  attributed  to  southward  propagation  of  the  rift. 

Chapman  &  Pollack  (1977)  suggest  that  heat-flow 
anomalies  in  Zambia  are  due  to  a  thinned  lithosphere  in  west 
and  central  Zambia.  Additionally,  incipient  rifting  may  be 
the  cause  of  a  long  wave  length  negative  Bouguer  gravity 
anomaly  to  the  SWS  of  Lake  Tanganyika;  similar  negative 
gravity  anomalies  that  coincide  with  surface  rift  features  in 
East  Africa  have  been  attributed  to  regions  of  thinning  in 
the  lithosphere. 

Maasha  &  Molnar  (1972)  have  shown  that  stress  drops 
outside  the  rift  in  southern  Africa  suggest  higher  stress 
there  than  inside  the  rift  to  the  north.  They  also  point  out 
that  volcanism  is  older  in  Ethiopia  than  in  Kenya,  suggesting 
that  the  Eastern  rift  may  have  grown  southward  during  its 


East  African  waveform  inversion  &  tectonics 

high  angle  of  faulting  for  this  event  shows  that  high 
deviatoric  tensional  stresses  continue  to  great  depths  and 
abnormally  high  crustal  temperatures  are  reflected  in  the 
region's  high  heat  flow.  It  is  very  likely  that  this  event 
occurred  below  the  brittle/ductile  transition  zone  which 
would  make  high  strain  rates  a  necessity. 

Meissner  &  Strehlau  (1982),  Strehlau  &  Meissner  U982), 
and  Strehlau  (1985)  show  that  very  high  strain  rates  are  a 
necessity  for  this  event.  Strehlau  &  Meissner  (1982)  suggest 
that  it  is  the  temperature  and  not  the  >  atrological 
properties  that  determines  the  cut-off  depth  of  seismicity; 
the  orientation  of  the  principle  stresses  and  the  subsequent 
fault  mechanism  also  play  a  minor  role  in  determining  the 
depth  of  seismic  and  aseismic  activity.  Chapman  &  Pollack’s 
(1977)  values  for  heat  flow  (q=7.1-  10~2  W/m)  =  1.69  HFU), 
thermal  gradient  (dT/dz-2. 54-  10-2  K/m) ,  and  thermal 
conductivity  (k=2.81  W/mK)  would  put  the  temperature  at  28  km 
above  600° C.  This  value  is  greater  then  the  critical 
seisraic/aseismic  transition  temperature  for  normal  faulting 
obtained  by  Strehlau  &  Meissner  —  500°C  q90° .  This  event, 
subsiquently ,  is  located  well  within  Strehlau  &  Meissner’s 
"plastic  domain"  where  stress  release  by  bulk  plastic  flow 
should  prevent  seismic  activity.  Based  on  these  simple 
models,  very  high  strain  rates  are  required  for  failure  at  23 
km. 

These  seismic  and  heat  flow  data  suggest  that  tensile 
stresses  extending  through  the  crust  are  separating  the  upper 
lithosphere  to  allow  the  upwelling  of  hot  mantle  material. 
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The  full  36  km  Gumper  &■  Pomeroy  crustal  model  was  use:? 
in  the  inversion  and  synthetic  computations  f^r  both  Zambia 
events  A  thicker  crust  than  that  used  for  the  Ethiopia  or 
Congo  events  is  consistent  with  plate  tectonic  assumptions 
for  this  area  because  well  defined  rifting  or  crustal 
thinning  is  not  immediately  evident.  Additionally,  the  lack 
of  well  defined  source  crustal  reverberation  phases  makes 
accurate  crustal  thickness  determination  for  either  event, 
rather  tenuous. 


South  Africa 


No  one,  as  of  yet,  has  suggested  that  the  rift-  extends 
the  whole  length  of  the  continent.  This  would  have  to  be  the 
case  if  the  South  African  event  of  September  29,  1969  is 
considered  to  be  due  to  southward  rift  migration  (figure  14). 
The  geology  of  the  Ceres  area  does  not  provide  substantial 
evidence  regarding  stress  orientation;  most  of  southern 
Africa  is  a  stable  shield  area  and  earthquakes  the  size  of 
the  September  29  event  (6.3)  are  rare.  In  fact,  no  previous 
earthquake  in  South  African  recorded  history  has  caused  such 
extensive  damage  (Green  &  Bloch  1971). 

An  oddity  of  this  event  is  that  the  fault  plane  from  the 
main  shock  obtained  by  Green  &,  McGarr  (1972)  does  not 
coincide  with  the  plane  delineated  by  the  carefully  studied 
aftershocks  (Green  &  Bloch  1971).  Our  corresponding  plane, 
however,  strikes  approximately  124® ,  only  5°  different  than 
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the  aftershock  data.  Green  &  McGarr  (1972)  report  that  there 
are  no  inconsistent  data  and  that  both  nodal  planes  are  well 
constrained. 

Inversion  yielded  a  nearly  pure  sinistral  strike  slip 
fault,  striking  at  about  124°  --  there  is  no  reported 
evidence  of  faulting  at  the  surface.  The  strike  of  this 
fault  coincides  nicely  with  that  of  the  Worcester  fault  zone 
which  apparently  has  not  been  active  since  the  Cretaceous 
(Green  &  McGarr  1972). 

The  P  waves  from  this  event,  as  previously  mentioned, 
could  not  be  consistently  modeled  using  simple  plane  layered 
structure.  The  gross  crustal  structure  in  the  area  of  the 
source  is  dominated  by  a  seaward  thinning  crust  to  the 
SSW-SE,  and/or  the  Karroo  basin  to  the  NE.  Figure  15 
illustrates  how  scattering  caused  by  a  dipping  interface 
distorts  P  waves  recorded  at  stations  on  opposite  sides  of  a 
node;  we  used  a  31  km,  3  layer  crustal  model  with  planar 
layers  and  a  seaward  rising  Moho  (10°).  While  this  figure 
does  not  explain  the  wave  form  distortions  recorded  at  the 
stations  to  the  NE,  it  does  show  how  dramatically  a  simple 
scatterer  affects  P  waves  radiated  from  a  strike  slip  source. 
The  curved  structure  of  the  Karroo  basin,  located  to  the  NE, 
may  be  the  cause  of  P  wave  distortions  observed  at  stations 
in  that  direction.  Since  there  are  few  studies  of  crustal 
heterogeneities  in  this  area  wc  cannot  determine  a  causal 
structure  which  produces  the  exact  waveform  changes  observed. 
However,  this  analysis  suggests  that  lateral  heterogeneities 
in  the  source  crust  are  a  plausible  explanation  of  the 
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effects  observed,  particularly  considering  the  inconsistent  P 
and  SH  mechanisms  found  in  the  inversions*.  It  is  also 
conceivable  that  multiple  fault  planes  were  involved  in  the 
rupture  causing  complex  waveforms. 


Summary  &  Conclusions 

All  of  the  mechanisms  obtained  in  this  study  ar£j 
somewhat  different  than  mechanisms  determined  by  previous 
workers  using  observational  and  P  wave  first  motion  data.  Our 
normal  fault  mechanisms,  for  the  Congo  and  both  Zambia 
events,  are  more  reliable  then  previous  interpretations 
because  of  the  inclusion  of  waveform  constraints.  The 
additional  information  provided  by  the  inversion,  focal  depth 
in  particular,  allows  for  greater  insight  into  the  tectonic 
activity  occuring  in  the  area  of  each  event. 

Eastern  Africa  has  been  shown  to  be  under  extensional 
stress  from  the  Afar  region  to  approximately  25°S.  The  fault 
plane  striking  327°  and  dipping  69° NE  has  been  chosen  as  the 
primary  node  for  the  Ethiopia  event  based  on  observational 
data.  This  mechanism  shows  that  the  Afar  area  is  under 
tensional  stress,  and  spreading  in  a  NNE-SSW  direction.  The 
Congo  mechanism  substantiates  the  accepted  view  that  this 
area  is  undergoing  crustal  extension.  Abundant  and  sundry 
geophysical  evidence  is  making  the  southward  rift  extension 
theory  nearly  inescapable.  The  large,  deep,  high  angle 
normal  fault  of  May  15,  1968  (Zambia)  shows  that  stress 
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accompanied  by  high  strain  rates  extends,  at  the  least, 
throughout  the  entire  crust.  Dist^  ;ted  P  waves  from  the 
nearly  pure  strike  slip  September  29,  1969  South  Africa  event 
demonstrate  the  existence  of  complex  crustal  structure  in  the 
area  of  the  event. 
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TABLE  1.2 


■time  functions  (liffp  )  -  1  second  box  cars 
Time  progresses  from  left  to  right,  top  to  bottom. 
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TABLE  1.3 


source  depth  and  condition  number 
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Ethiopia 
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30 

Congo 

10 

82 
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28 
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28 

26 

figure  2  -  Lower  hemisphere  projections  of  mechanisms  determined  by 
previous  authors  and  this  study;  P=pressure  axis,  T=tension  axis. 

1:  Ethiopia  event  -  (a)  McKenzie  et  al.  (1970),  (b)  this  study.  2: 
Congo  event  -  (a)  Sykes  (1967),  (b)  this  study.  3:  Zambia  (December) 
-  (a)  Fairhead  &  Girdler  (1971),  (b)  this  study.  4:  Zambia  (May)  - 
(a)  Fairhead  &  Girdler  (1971),  (b)  this  study.  5:  South  Africa  -  (a) 
Maasha  &  Molnar  (1972),  (b)  this  study,  (c)  Fairhead  &  Girdler  (1971), 
(d)  Green  &  McGarr  [shows  range  of  fault  planes]  (1972). 
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figure  5  -  Lower  hemisphere  projections  [(a)  P  waves,  (b)  SH  waves]  of 
fault  mechanism  for  the  Congo  event  of  March  20,  1966;  +  signs  represen 
compression,  and  -  dilatation.  The  pressure,  tension  and  null  axes  are 
represented  by'P,  T  and  1  respectively.  The  thicker,  darker  traces  are 
those  of  the  observed  dat* . 
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figure  7  -  Lower  hemisphere  projections  [(a)  P  waves,  (b)  SH  waves]  of 
fault  mechanism  for  the  Zambia  event  of  May  15,  1968;  +  signs  represent 
compression,  and  -  dilatation.  The  pressure,  tension  and  null  axes  are 
represented  by  P,  T  and  I  respectively.  The  thicker,  darker  traces  are 
those  of  the  observed  data. 
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figure  8a 

figure  8  -  Lower  hemisphere  projections  [(a)  P  waves,  (b)  SH  waves]  of 
fault  mechanism  for  the  Zambia  event  of  December  2,  1968;  +  signs  represent 
compression,  and  -  dilatation.  The  pressure,  tension  and  null  axes  are 
represented  by  P,  T  and  I  respectively.  The  thicker,  darker  traces  are 
those  of  the  observed  data. 
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figure  10  -  Histogram  representing  the  number  of  short  period  records 
with  a  given  sP-P  time  for  he  September  29,  1969  South  Africa  event 
The  short-period  P  wave  ro  rded  at  ARE  (WWSSN  station)  shows  sP-P 
separation.  Measurement  error  is  estimated  to  be  qO.l  seconds. 
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(b) 

figure  13  -  (a)  -  Epicenter  location  (black  star)  and  faults  in  the  area  of  the 
March  20,  1966  Congo  event.  Hatched  area  in  the  upper  right  is  Lake  Albert,  and 
in  the  lower  left  is  Lake  Edward  (Lahr  &  Pomeroy,  1970).  This  area  is  Contained 
within  block  2a  of  figure  1.  (b)  -  Epicentral  locations  of  the  Zambia  events, 
map  from  Chapman  &  Pollack  (1977).  Stars  represent  the  location  of  the  December 

2nd  (upper  left),  and  May  15th  (lower  right)  events.  This  area  is  contained 
within  block  2b  of  figure  1. 
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figure  14  -  Faulting  in  the  area  of  the  September  29,  1969  South  African 
event;  the  large  (M=5.9)  dot  represents  the  main  shock  (Fairhead  & 
Girdler,  1971).  This  area  is  contained  within  block  3  of  figure  1. 
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figure  15  -  Observed  P  wave  data  from  the  September  29,  1969  South 
Africa  event  plotted  adjacent  to  synthetics  computed  assuming  a 
seaward  thinning  crust  (dipping  moho) . 


SECTION  7 


The  Pennsylvania  State  University 
The  Graduate  School 
Department  of  Geosciences 


Upper  Mantle  P-Velocity  Structure 
between  Eastern  North  America  and  Hispaniola 


A  Paper  in 
Geophys i cs 

by 

Nancy  L.  Niemann 


Submitted  in  Partial  Fulfillment 
of  the  Requirements 
for  the  Degree  of 


Master  of  Science 
December  1986 

©  1986  by  Nancy  L.  Niemann 


I  grant  The  Pennsylvania  State  Uni  vers  it.'  the 
nonexclusive  right  to  use  this  work  for  the  University 
own  purposes  and  to  make  single  copies  of  the  work 
available  to  the  public  on  a  not-for-profit  basis  if 
copies  are  not  otherwise  available. 


Nancy  L.  Niemann 


ACKNOWLEDGMENTS 


I  would  like  to  express  my  gratitude  to  and 
admiration  for  Dr.  Charles  A.  Langston,  who  proposed 
this  study,  provided  the  software  that  calculated  the 
WKBJ  seismogram,  and  advised  me  throughout  my  research. 
His  insight,  experience,  and  above  all  patience  allowed 
the  completion  of  this  work. 

To  David  Johnston  go  my  warm  thanks  for  making 
available  his  digitized  records  of  the  Hispaniola  earth¬ 
quake.  Richard  "Boomer"  Baumstark  at  the  Center  for 
Seismic  Studies,  Washington  DC,  enabled  my  access  to  the 
RSCP  data  set. 

Highest  praises  go  to  the  Geosciences  Department , 
and  to  Dr.  Langston  in  particular,  for  procuring  and 
maintaining  the  PRIME  minicomputer,  on  which  I  made  my 
calculations.  I  used  the  word  processing  capabilities 
of  the  Geophysics  Program  and  of  Conoco,  Inc.,  to 
prepare  this  manuscript. 

My  research  was  partially  supported  by  the  Advanced 
Research  Projects  Agency  and  monitored  by  the  Air  Force 
Office  of  Scientific  Research  under  Contract 
F49620-83-K-0019.  The  Geosciences  Department  provided 
additional  support  with  assistantships  for  teaching  and 
research . 


INTRODUCTION 


Seismic  waveforms  recorded  at  stations  15°  to  30° 
from  earthquake  sources  have  produced  as  much  knowledge 
about  the  methods  used  to  study  them  as  they  have  about 

the  P  and  S  velocities  of  the  upper  mantle.  In  the 

present  study,  an  attempt  to  define  a  slight  lateral 
variation  in  the  depth  of  the  650-km  discontinuity 
contributes  to  this  knowledge. 

The  eastern  Hispaniola  earthquake  of  September  14, 
1981,  with  body-wave  magnitude  5.9  (International 
Sei smol ogical  Centre,  1983),  generated  waveforms  that 
are  useful  for  studying  wave  propagation  beneath  the 
ocean  off  the  southeast  coast  of  the  United  States.  The 

data  used  in  this  study  are  adequate  to  study  the  region 

defined  by  +'s  in  Figure  1.  The  earthquake  source 
function  is  simpler  than  many  that  have  been  modeled  in 
previous  upper  mantle  waveform  studies  (cf.  Burdick, 
1981).  The  event's  hypocenter  was  172+1.9  km  deep 
(International  Seismological  Centre,  1983).  Phases 
arriving  after  direct  P  are  sufficiently  delayed  that 
they  do  not  interfere  with  the  triplicated  P  phases.  In 
principle,  therefore,  the  simpler  source  function  and 
greater  depth  should  make  observing  the  triplication 
branches  easier. 
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We  attempt  to  interpret  the  P  waveforms  in  the 
first  5  to  10  seconds  of  the  upper  mantle  records  by 
comparing  them  to  synthetic  seismograms  calculated  for 
Burdick's  (1981)  S8  model  for  an  adjacent  and  partially 
overlapping  region.  Since  S8  does  not  predict  the 
Hispaniola  travel  times  and  waveforms  well,  we  attempt 
to  find  a  new  model  that  does. 

Johnson  (1967)  inverted  dT/dA  data  to  obtain  a 
P-velocity  model  for  the  western  North  American  tectonic 
region.  Helmberger  and  Wiggins  (1971)  began  their  upper 
mantle  studies  by  inverting  travel-time  data  to  derive  a 
starting  velocity  model.  They  proceeded  to  manipulate 
the  model  by  trial  and  error  until  waveforms  calculated 
from  it  resembled  observed  waveforms. 

If  a  published  velocity  model  exists,  then  invert¬ 
ing  travel  times  to  obtain  a  starting  model  is  not  nec¬ 
essary.  A  number  of  authors  have  successfully  manipu¬ 
lated  published  starting  models  derived  by  waveform 
analysis  for  regions  that  they  were  studying  with  more 
recently  available  data  (Helmberger,  1973,  Wiggins  and 
Helmberger,  1973,  Helmberger  and  Engen,  1974,  Burdick 
and  Helmberger,  1978,  Given  and  Helmberger,  1980, 

Burdick,  1981,  and  Grand  and  Helmberger,  1984).  The 
availability  of  model  S8  obviates  inverting  travel-time 
data  to  obtain  a  starting  velocity  model  here. 

Burdick  (1981)  set  out  to  demonstrate  that  varia¬ 
tions  in  the  upper  mantle  P-velocity  structure  beneath 


STUDY  AREA 


The  Hispaniola  earthquake  of  September  14,  1981, 
was  recorded  at  stations  in  Bermuda  and  in  eastern  North 
America  at  ranges  between  14°  and  28°  (Figure  1).  The 
region  of  study  is  defined  by  the  distribution  of  the 
bottoming  points  ( + '  s  in  Figure  1}  of  the  source-to- 
station  ray  paths.  The  range  in  source-to-station  azi¬ 
muths  of  these  stations  is  abort  65°.  Naturally  the 
character  of  the  wave  propagation  paths  from  the  source 
to  each  station  will  change  with  azimuth  and  with 
di stance . 

Upper  mantle  velocity  studies  are  traditionally 
conducted  for  geographic  regions  for  which  large  numbers 
of  upper  mantle  observations  have  accumulated.  Clearly, 
the  least  ambi guous  solution  to  the  inversion  problem 
will  result  if  the  observations  are  as  evenly  and 
closely  distributed  with  distance  as  possible.  We  must 
balance  the  advantages  of  working  with  a  large  data  set 
against  the  disadvantage  of  averaging  too  large  a 
region,  since  the  larger  the  area,  the  greater  the 
chance  of  including  regions  of  the  upper  mantle  with 
slightly  varying  charactei. 


Previous  authors  (Helmberger  and  Wiggins,  1971, 
Wiggins  and  Helinberger,  1973  ,  and  Burdick,  1981)  have 
found  that  describing  the  upper  mantle  data  of  North 
America  requires  two  velocity  models:  one  for  the 
tectonic  area  to  the  west  and  a  second  for  the  stable 
shield  region.  If  the  upper  mantle  beneath  the  stable 
continental  shield  region  has  little  lateral  variation, 
so  that  source-to-station  wave  propagation  paths  are 
similar,  then  one  model  should  fit  most  of  the  data. 

All  of  the  North  American  stations  at  upper  mantle 
distances  from  the  Hispaniola  event  recorded  waves  that 
had  traveled  paths  underneath  first  oceanic  crust  and 
then  continental  shelf  and  stable  continental  shield 
crust;  each  of  the  travel  paths  bottoms  beneath  oceanic 
crust.  Hodel  S8  was  developed  fro'.i  source-to-station 
ray  paths  that,  in  almost  every  case,  sampled  the  upper 
mantle  beneath  stable  continental  areas.  North  American 
records  of  the  1978  Bermuda  earthquake  are  a  subset  of 
this  data  set.  Like  Hispaniola  source-to-station  paths, 
those  from  the  Bermuda  event  are  partly  oceanic  and 
bottom  beneath  oceanic  crust.  Since  the  region  covered 
by  the  Hispaniola  stations  overlaps  the  Bermuda  region 
from  the  south,  upper  mantle  wave  propagation  should  be 
similar  in  both  cases.  Since  the  S8  model  averages 
European  and  North  American  upper  mantle  data,  however, 
it  may  not  fit  the  Hispaniola  data  very  well. 


Reduced  Travel  Time  (sec) 


Figure  3a.  Reduced  travel-time  curve  for  the  S8  upper 
mantle  velocity  model.  Travel  times  are  reduced  by 
10.3  km/sec. 
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Figure  3b.  Reduced  travel-time  curve  for  the  E19  upper 
mantle  velocity  model.  Travel  times  are  reduced  by 
10.3  km/sec. 
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samples  signals  at  10  sampl es/second ,  proved  to  be  the 
single  most  useful  data  record  available  to  this  study. 
With  proper  adjustments  to  the  S8  model,  the  synthetic 
seismogram  that  we  calculated  for  SCP  closely  resembles 
the  intermediate-period  SCP  data  in  waveform  shape  and 
smoothness . 


SYNTHETIC  SEISMOGRAM  METHODS 


It  is  reasonable  to  consider  a  far-field  seismogram 
as  a  convolution  of  a  number  of  operators.  For  example, 
the  vertical  displacement,  R(t),  is  given  by: 

R(t)  =  E ( t )  *  I ( t)  *  A ( t )  *  S(t)  *  M(t), 
where  E(t)  is  the  earthquake  source  time  function,  I ( t ) 
is  the  instrument  response,  A ( t )  is  the  attenuation 
operator,  S(t)  represents  the  local  source  and  receiver 
effects,  M(t)  is  the  mantle  response  operator,  and  * 
denotes  convolution. 

The  mantle  response  operator  M(t)  consists  of  the 
set  of  Green's  functions  for  waveforms  propagating  from 
the  Hispaniola  event  through  the  mantle  velocity  models. 
Johnston  and  Langston  (1984)  determined  the  earthquake 
source  parameters  for  the  Hispaniola  event  by  moment 
tensor  inversion;  M(t)  includes  the  fault-plane  orienta¬ 
tion,  seismic  moment,  and  source  depth  of  Model  C  of 
Johnston  and  Langston  (1984).  The  WKBJ  seismogram 
method  of  Chapman  (1978)  is  used  to  calculate  the  mantle 
response.  Advantages  and  disadvantages  of  calculating 
upper  mantle  seismograms  with  WKBJ  theory  are  documented 
by  Grand  and  Helmberger  (1984)  and  Given  (1984). 


This  method  assumes,  of  course,  that  at  teleseismic 
distances  the  earth's  response  to  a  propagatingwave  is 
a  delta  function  which  has  been  smoothed  by  attenuation. 
We  also  assume  that  the  teleseismic  waveforms  are  not 
affected  by  receiver  complications.  The  observation 
that  the  Hispaniola  teleseismic  vert i ca 1 -componen t  wave¬ 
forms  (see  Figure  4)  are  remarkably  similar  at  different 
azimuths  and  distances  validates  these  assumptions. 
Moreover,  Given  (1984)  and  Grand  and  Helmberger  (1984) 
demonstrate  that  convolving  a  representative  teleseismic 
wavelet  with  mantle  impulse  responses  gives  upper  mantle 
synthetics  that  are  comparable  to  waveforms  calculated 
using  theoretical  operators. 

Several  long-period  teleseisms  for  the  Hispaniola 
event  are  shown  in  Figure  4.  During  the  first  20 
seconds  of  each  record,  two  distinct  pulses  arrive. 
Direct  P  is  consistently  followed  after  about  7  seconds 
by  an  event  with  reversed  polarity.  Johnston  and 
Langston  (1984)  interpreted  the  second  event  as  an  up- 
going  ray  from  the  source  that  reflects  off  the  bottom 
of  the  low-velocity  zone,  then  travels  a  path  similar  to 
that  of  P.  To  make  this  interpretation,  they  recon¬ 
structed  the  waveforms  by  convolving  the  earth  response 
to  the  two  rays  (calculated  using  generalized  ray 
theory)  with  the  appropriate  instrument,  earthquake 
source  time  function,  and  attenuation  operators.  The 
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3  seconds  long,  as  the  source  mechanism  study  of 
Johnston  and  Langston  (1984)  determined. 


RECEIVER  FUNCTION 


In  this  study,  as  in  other  upper  mantle  waveform 
studies,  a  delta  function  represents  the  effect  of 
shallow  structure  beneath  the  receiver.  The  receiver 
structure  is  generally  not  known.  Assuming  that  it  has 
no  effect  on  wave  propagation  in  this  case  simplifies 
the  effort  of  calculating  synthetic  seismograms. 

Of  course,  this  assumption  is  valid  only  for 
stations  whose  crustal  structure  does  not  create  arri¬ 
vals  that  are  obvious  in  the  P  waveforms.  It  requires 
stations  at  which  converted  phases  and  reverberations  in 
the  crust  are  not  strong  enough  to  complicate  the  upper 
mantle  waveforms.  Of  the  WWSSN  stations  that  recorded 
the  Hispaniola  event,  Burdick  and  Langston  (1977)  sug¬ 
gest  that  only  MNT  (Montreal,  Canada)  may  be  underlain 
by  crust  with  sharp  velocity  contrasts  that  cause  impor¬ 
tant  crustal  phases. 

To  determine  whether  a  delta  function  closely 
approximates  the  receiver  function  of  each  of  the  upper 
mantle  stations  used  here,  the  reasoning  of  Burdick  and 
Langston  (1977)  applies.  The  presence  of  SV  energy, 
indicating  that  crustal  phases  are  interfering  with  the 
vertical  P  waveform,  is  evident  on  the  rad i a  1 -component 
record.  If  the  P  waveforms  on  the  vertical-  and  radial- 


component  records  are  similar,  then  the  receiver  crjst 
is  "transparent"  {Burdick  and  Helmberger,  1978,  p  1703) 
to  propagating  waves,  and  the  single  delta  function  is 
appropriate.  If  the  horizontal  components  do  not  resem¬ 
ble  the  vertical,  then  the  crustal  structure  contains 
sharp  velocity  contrasts  whose  effects  on  wave  propa¬ 
gation  cannot  be  ignored.  In  such  a  case,  the  compli¬ 
cations  in  the  vertical  P  waveform  are  not  due  simply  to 
upper  mantle  structure. 

RSCP's  horizontal  components  were  rotated  into 
radial  and  transverse  directions.  MNT,  OTT,  and  PAL  are 
all  within  11°  of  due  north  of  the  source;  thus,  they 
are  nearly  naturally  rotated.  Radial  and  vertical  com¬ 
ponents  in  each  case  do  not  differ  markedly.  The  trans¬ 
verse  components  are  of  very  low  amplitude,  indicating 
that  little  SH  energy  is  generated  in  the  receiver 
crust. 

For  the  purpose  of  this  study,  digitizing  and  rota¬ 
ting  the  horizontal  components  for  the  remaining  sta¬ 
tions  were  unnecessary.  The  NS  and  EW  components  for 
each  station  (except  SCP)  are  similar  enough  to  each 
other,  and  to  the  vertical  component,  that  no  station 
was  discarded. 

The  NS-component  intermediate-period  record  for  SCP 
is  lot  available.  While  the  long-period  data  encourage 
us  that  SCP's  crustal  structure  is  acceptable,  crustal 
phases  should  be  less  obvious  on  long-period  records 
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than  on  short-period  records  (Burdick  and  Kelmberger, 
1978,  and  Given  and ’ Helmberger ,  1980).  For  the  first  3 
to  5  seconds,  the  intermediate-period  EW  record  for  SCP 
is  identical  to  the  vertical,  however,  indicating 
insignificant  crustal  features.  In  support  of  this, 
Langston  and  Isaacs  (1981)  successfully  modeled  the 
crust  and  uppermost  mantle  beneath  central  Pennsylvania 
as  a  41-km  thick  homogeneous  layer  over  a  halfspace. 
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TECHNIQUE  AND  RESULTS 


Analysis  of  upper  mantle  P  waveforms  is  at  the  heart 
of  this  study..  We  interpret  an  observed  waveform  by 
constructing  a  synthetic  meant  to  predict  features  in 
the  waveform  and  then  by  comparing  the  fit  between  the 
two.  Determining  the  quality  of  the  fit  of  predictions 
to  observations  is  partly  subjective  in  this  study. 

This  means  that  we  must  define  the  characteristics  of  a 
synthetic  waveform  that  is  a  "good  fit"  to  an  observa¬ 
tion.  Of  course  the  synthetics  are  only  approximations 
to  the  observations  and  are  limited  by  the  assumptions 
that  we  have  made.  For  example,  we  have  attenuated  the 
long-period  synthetics  too  much  and  so  have  removed  the 
high  frequencies  that  are  present  in  the  observations. 

We  pay  attention  then  to  matching  the  features  of  the 
waveforms  that  are  evidence  of  triplicated  direct  P 
arrivals:  the  travel  time  of  the  first  arrival,  the 

relative  timing  of  subsequent  arrivals,  the  relative 
amplitudes  of  the  pulses,  the  widths  of  the  pulses,  and 
the  polarities  of  the  pulses.  Quantitatively,  relative 
travel  times  that  are  within  one  second  of  those  in  the 
observations  and  relative  amplitude  measurements  that 
are  within  152  of  those  in  the  observations  are  accep¬ 
table. 
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Figure  9.  SCP  i  ntermedi ate- peri od  data  an< 
Pulses  represent  arrivals  of  rays  that  1)  ' 
below  the  670-km  discontinuity,  2)  turned  , 
and  3)  traveled  through  it.  The  starred  ai 
data  is  not  modeled  by  our  assumptions. 
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Figure  11.  In  termed  i  a  te-per-1  od  synthetic  seismograms 
calculated  for  S8.  Absolute  travel  times  are  reduced  by 
10.3  km/sec  and  advanced  by  30  seconds.  Note  that  the 
SCP  observation  (figure  9)  resembles  the  S8  synthetic  at 
22.25° 


reason  we  attempt  to  modify  the  lower  half  of  the  velo¬ 
city  profile.  We  suggest  that  by  shifting  the  6 5 0 - km 
triplication  (DEF  in  Figure  3a)  along  the  10.3-km/sec 
branch,  we  will  not  only  systematically  delay  travel 
times,  but  we  will  also  change  the  relative  timing  of 
arrivals. 

At  a  point  to  the  right  of  F  in  Figure  3a,  the  shift 
will:  1)  delay  all  arrivals,  2)  change  the  relative 

timing  of  arrivals  and,  in  particular,  advance  the  CD 
branch  with  respect  to  the  EF  branch,  and  3)  compress 
the  entire  triplicated  pulse. 

To  shift  and  reshape  the  650-km  triplication 
branches,  we  will  adjust  the  starting  velocity  profile 
in  the  following  ways: 

1)  lower  the  650-km  discontinuity  to  shift  the 
triplication  out  in  distance, 

2)  vary  the  gradient  above  and  below  the 
discontinuity  and  change  the  absolute  size  of 
the  discontinuity  to  distort  the  triplication, 
and 

3)  sharpen  or  smooth  the  discontinuity  to  lengthen 
or  shorten  the  duration  of  the  triplication. 

We  manipulate  the  starting  model  according  to  these 
guidelines  until  we  achieve  a  velocity  profile  whose 
synthetics  satisfactorily  match  the  data.  The  resulting 
model  will  represent  only  one  combination  of  these  acti¬ 
vities,  however.  Since  other  combinations  will  shape 
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profiles  that  predict  the  observations  as  weii,  we  con¬ 
clude  that  our  model  is  not  unique. 

Our  strategy  is  to  adjust  the  model  so  that  it  fits 
the  SCP  intermediate-period  (IP)  record  (Figure  9).  At 
this  range,  direct  P  arrivals  are  separated  by  about  1 
to  2  seconds.  A  time  separation  of  this  order  is  resol- 
val  on  tho  'P  record.  This  record  is  so  filtan  and 
noise-free  that  we  can  calculate  synthetics  that  almost 
duplicate  it.  Once  we  have  established  this-  one  data 
point  on  the  travel-time  curve,  we  proceed  to  fit  the 
remaining  observations  as  closely  as  possible.  Represen¬ 
tative  waveforms  for  the  new  model  are  shown  in  Figure  6. 
S8  and  the  new  velocity  model  E 1 9  are  listed  in  Table  1 
and  plotted  in  Figure  2. 

The  400-km  Discontinuity.  Stations  located  between  12° 
and  20°  from  the  source  are  in  an  excellent  position  to 
record  triplicated  P  arrivals  from  the  vicinity  of  the 
400-km  discontinuity.  Of  the  stations  whose  records  are 
available  for  this  study,  only  BEC  at  14.5°  falls  in 
this  range.  The  S8  model  predicts  a  waveform  that  fits 
the  BEC  observation  acceptably  well.  The  lack  of  obser¬ 
vations  from  this  distance  range  is  troubling,  however. 
Fortunately,  observations  of  the  AB  branch  (Figures  3a 
and  3b)  beyond  21°  help  to  confirm  that  the  triplication 
caused  by  the  upper  discontinuity  in  the  S8  velocity 
model  satisfies  the  data. 
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Table  1.  S8  and  E19  P  Velocity  Models. 


Depth 

S8 

E 19 

30.0  km 

6.50  km/sec 

6.500  km/sec 

50.0 

7.10 

7  . 100 

60.0 

8.10 

8.100 

65.0 

8.20 

8.200 

75.0 

8.24 

8.240 

135.0 

8.26 

8.260 

175.0 

8.30 

8.300 

245.0 

8.44 

8.440 

390.5 

8.81 

8.810 

400.0 

9.25 

9.250 

600.0 

10.000 

615.0 

10.07 

625.0 

10.11 

10.175 

645.0 

10.25 

650.0 

10.76 

660.0 

30.225 

667.5 

10.275 

672.5 

10.350 

675.0 

10.400 

675.1 

10.650 

680.0 

10.700 

690.0 

10.775 

705.0 

10.850 

805.0 

11.060 

905.0 

11.26 

11.260 

V' VV 


TOjvtV 


wVvl  V VsA1 OlVC%  '.‘"v  s  's.YxV* 


From  21°  to  24°,  the  spread  between  the  AB  and  CD 
branches  of  the  triplication  increases  from  about  5 
seconds  to  about  10  seconds.  The  position  of  the  AB 
branch  in  this  range  represents  the  travel  time  of  a 
diffracted  arrival  from  just  above  the  400-km  disconti¬ 
nuity.  The  WKBJ  method  predicts  too  large  an  amplitude 
for  this  pulse,  which  arrives  beyond  the  B  cusp  of  the 
travel-time  curve.  Despite  this,  the  pulse's  travel 
time  is  useful  information.  Before  the  pulse  from  the 
400-km  discontinuity  arrives,  of  course,  we  see  arrivals 
from  the  deeper  discontinuity.  At  BLA  and  SHA  the  arri¬ 
val  from  the  region  of  the  400-km  discontinuity  causes 
the  second  peak  in  the  waveform.  By  23°  at  PAL  and  23.7° 
at  SCP  (long-period  record),  the  amplitude  of  the 
arrival  has  decreased  abruptly  and  the  AB  branch  has 
moved  back  in  time.  The  disappearance  of  the  AB  branch 
is  evidence  that  the  first  triplication  should  be  no 
1 arger  than  it  is. 

The  6 70- km  Discontinuity.  Our  model  predicts  observed 
absolute  travel  times  to  within  4  seconds  at  all  sta¬ 
tions.  We  overestimate  the  travel  time  at  BEC  by  2 
seconds,  at  BLA  by  2.5  seconds,  at  SCP,  AAM,  FVM,  MNT, 
and  OTT  by  4  seconds,  and  underestimate  at  RSCP  by  2 
seconds.  The  SHA  and  PAL  first  arrivals  are  on  time; 
perhaps  the  blame  for  inconsistency  of  the  SHA  time, 
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which  should  be  closer  to  the  b L A  time,  lies  in  error  in 
digitizing  the  record,  as  we  mentioned  earlier. 

Waveform  fits  to  BLA  and  SHA  are  very  good.  The 
first  cycle  of  the  waveform  is  slightly  broader  than  the 
source  wavelet.  The  pulse  from  just  above  the  670-km 
discontinuity  is  followed  after  one  second  by  the  pulses 
that  propagate  through  the  discontinuity  and  along  the 
underside  of  it  (Cormier  and  Choy,  1981).  The  second 
peak  is  matched  closely  by  a  pulse  that  we  predict 
arrives  from  the  vicinity  of  the  400-km  discontinuity. 

Our  assumptions  are  not  complete  enough  to  explain  the 
slight  discrepancy  in  relative  amplitudes  between  SHA 
and  BLA.  The  discrepancy  may  be  due  to  lateral  variation 
in  mantle  structure.  T.ie  azimuthal  separation  between 
these  two  stations  with  respect  to  the  source  is  25°. 

Across  the  CD-EF  intersection  at  SCP,  the  interme¬ 
diate-period  synthetic  fits  the  observation  almost  exact¬ 
ly.  The  travel-time  curve  indicates  a  time  spread  be¬ 
tween  the  EF  branch  and  the  back-branch  of  slightly  less 
than  2  seconds;  this  is  precisely  the  time  spread  between 
the  first  trough  and  the  second  peak  in  the  data.  The 
onset  of  the  pulse  from  the  CD  branch,  0.8  seconds  behind 
the  EF  branch,  is  distinct  in  the  short-period  digital 
record  (Figure  12).  This  confirms  the  time  separation 
between  the  first  two  arrivals.  More  importantly,  the 
fact  that  no  arrivals  intervene  between  the  two  supports 
our  interpretation  that  these  pulses  are  indeed 
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Figure  12.  Short-period  record  of  the  Hispaniola  event 
recorded  at  SCP  (DWWSSN).  The  first  pulse  (trough  at  1) 
arrives  from  below  the  670-km  discontinuity.  The  trough 
at  2  represents  a  ray  that  turned  above  the  discontinuity. 
The  peak  at  3  corresponds  to  an  arrival  on  the  back  branch 
(DE  in  figure  3b),  a  ray  that  turned  within  the 
discontinuity.  Relative  arrival  times  predicted  by  E19 
are  close  to  the  timing  of  arrivals  in  this  record. 


triplicated  direct  P  arrivals.  We  can  see  each  of  the 
three  triplication  branches  despite  slight  i nterference 
among  the  pulses.  The  two  pulses  with  negative  polarity 
are  arrivals  from  below  the  670-km  discontinuity  and 
from  just  above  it,  in  that  order.  Smoothing  the  corners 
of  the  discontinuity  in  the  profile  gave  the  pulses  the 
correct  relative  amplitudes.  The  strong  peak  is  the 
arrival  of  a  ray  that  traveled  through  the  discontinuity; 
this  is  the  only  record  in  which  the  90-degree  phase- 
shifted  pulse  from  the  back-branch  is  directly  observ¬ 
able.  The  hump  on  the  right  side  of  the  peak  ( *  in 
Figure  9)  is  not  modeled  by  our  assumptions. 

The  quality  of  waveform  fits  at  AAM,  FVM,  MNT,  and 
OTT  is  very  difficult  to  judge.  MNT  and  OTT  especially 
are  very  close  to  the  P  nodal  plane;  we  expect  the  trip¬ 
licated  P  arrivals  to  be  of  very  low  amplitude.  Never¬ 
theless,  they  are  distinguishable  from  noise,  and  vital 
in  confirming  the  position  of  the  CD  branch.  The  first 
trough  in  the  synthetics  corresponds  to  the  first  trough 
in  each  of  the  observations.  The  relative  timing  of 
peaks  and  troughs  that  our  model  predicts  agrees  well 
with  the  data.  The  relative  amplitudes  of  the  first 
trough  and  second  peak  are  almost  exact  in  each  case. 

The  amplitude  of  the  first  peak,  however,  is  slightly 
high,  especially  at  AAM  and  FVM.  The  first  trough  repre¬ 
sents  the  head  wave  that  travels  urderneath  the  670-km 
discontinuity.  The  second  pulse,  the  superposition  of 
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Figure  13.  RSCP  intermediate-period  data  and  synthetics. 
Neither  model  predicts  the  data  well.  The  synthetics 
appear  to  be  extremely  sensitive  to  the  position  of  the 
intersection  of  the  CD  and  EF  branches  of  the  travel-time 
curve.  We  had  no  success  in  modeling  this  observation. 
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International  Se i smologi cal  Centre,  1983)  would  bring 
predicted  times  closer  to  the  data.  This  shift  in  travel 
times  would  be  accomplished  without  changing  the  tripli¬ 
cations.  On  the  other  hand,  thinning  the  crust  in  the 
model  is  a  more  likely  solution.  The  crust  in  E19  is  60 
km  thick,  which  could  be  characteristic  of  some  stable 
continental  shield  regions  on  other  continents,  but  is 
too  thick  for  the  crust  beneath  SCP  and  likely  for  the 
other  stations  as  well.  Substituting  a  40-km  thick 
crustal  model  with  velocities  of  6.0  km/sec  for  the 
upper  20  km  and  6.8  km/sec  for  the  lower  20  km  would 
shorten  the  travel  time  through  the  crust  by  1.1 
seconds . 


r  <- .V ,  /:  y.  A  y.v. W.  A  v.  A  A  rfVX .  'Oiru.T-l* - 


306 


DISCUSSION 

Burdick's  (1981)  model  S8  has  been  modified  to 
satisfy  a  new  set  of  upper  mantle  data  from  eastern 
North  America.  The  new  model,  E19,  is  identical  to  S8 
above  about  640  km.  Features  that  distinguish  E19  from 
S8  are: 

1)  the  deepest  velocity  jump  occurs  at  670  km 
instead  of  at  S8's  650  km; 

2)  the  velocity  jump  is  not  discontinuous  as  in  S 8 ; 
the  velocity  jump  occurs  over  a  depth  range  of 
about  30  km. 

The  increased  amount  of  structure  resolvable  near 
the  670-km  discontinuity,  feature  2  above,  is  a  direct 
benefit  of  the  use  of  the  intermediate-period  SCP 
record.  This  record  was  sufficiently  noise-free  that  we 
were  able  to  reproduce  its  first  3  to  4  seconds  excep¬ 
tionally  well.  Greater  availability  of  such  high 
quality  broad-band  digital  data  would  no  doubt  enhance 
similar  studies. 

The  claim  of  Burdick  and  other  authors  that  mantle 
velocity  is  laterally  homogeneous  below  about  250  km 
raises  questions  about  the  validity  of  the  new  model. 
Does  disagreement  about  the  position  of  the  lower 
discontinuity  imply  inhomogeneity?  Is  a  difference  in 
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in  a  starting  velocity  model.  Wi 
attempt  to  answer  the  above  quest 


to  Given's  cannot  distinguish  differences  in  20  km  in 
the  position  of  a  change  in  gradient. 

The  fact  remains  that  the  new  model  fits  the 
Hispaniola  data  better  than  S8  does.  The  new  model  is 
not  unique,  however.  Given  (1984,  p.  55)  described  the 
uniqueness  problem  in  trial-and-error  inversions  in  this 
way:  "The  difficulty  of  finding  a  velocity  model  that 

satisfied  the  data  is  implicitly  used  as  an  argument 
that  t  ,ie  model  was  unique."  We  noted  that  particular 
combinations  of  gradient  steepness,  discontinuity  depth, 
and  velocity  jump  would  satisfy  the  data.  The  new  model 
is  one  combination.  Other  combinations  would  fit  the 
data  as  well.  Adding  more  observations  would  permit 
others,  until  we  accumulate  enough  data  so  that  we  may 
begin  to  better  constrain  the  model.  Even  though  the 
average  distance  spacing  of  Burdick's  observations  is 
0.5°  to  1.0°,  S8  is  not  unique  either. 
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